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PREFACE 


.This volume is Che first of a three volume set presenting 
the description and program documentation of a mathematical model 
package for thermal pollution analyses and prediction. This vol- 
ume presents Che mathematical formulation of these models, including 
assumptions, approximations, governing equations, boundary and in- 
itial conditions, numerical method of solution and sample results. 
The two model formulations are the rigid-lid and free-surface, re- 
spectively. These programs were developed by the Thermal Pollution 
Group at the University of Miami, and were funded by the National 
Aeronautics and Space Administration (NASA) , thus the program 
names NASUM I, NASUI4 II and NASUM III were given to reflect this 
joint effort. 

These models are three-dimensional and time-dependent using 
the primitive equation approach. They have sufficient generality 
in programming procedure to allow application at sites with diverse 
topographical, features. NASUM I is a rigid-lid formulation and is 
presented in detail in Volume II. NASUM I consists of both near 
and far field versions. The near field simulates themal plume 
areas, and the far field version simulates larger receiving aquatic 
ecosystems. The models in NASUM I simulate the velocity and tem- 
perature fields for given meteorological and plant intake and dis- 
charge conditions. Three versions of the rigid-lid formulation are 
presented in Volume II comprising NASUM I; one for near field 
simulation, the second for far field unstratified situations, and 
the third is for stratified basins, far field simulation. 


iTASUil II and NASUl-I III are Irec-ijurface lormulaLiona 
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and are presented in detail in Volume III, Both pro>>rams 
present surface height variations, velocity field and tem- 
perature field for the '’complete field", MASUl-I II is a 
far field formulation and ?'s used without including the plant 
thermal discharge. NASUM III used horizontal stretching in 
order to provide higher resolution at thejrmal discharge point, 
as well as including far field influences such as varying tide 
and ambient currents at points far from the point of discharge. 
It also includes far field influences such as varying tide 
and ambient currents at points far from the point of discharge. 

The three volumes are intended as User's manuals and, 
accordingly, they present specific instructions regarding data 
preparation for program execution and specific sample problems. 
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LIST OF SYMBOLS FOR RIGJD-LID MOnEL 


Aji Horizontal kinematic eddy viscosity 
Vertical kinematic eddy viscosity 
A, Vertical eddy viscosity 

IM 

^ref kinematic eddy viscosity 

Av 

Bj, Horizontal eddy thermal dlffuslvlty 

Vertical eddy thermal diffusivity 

®ref eddy ther:"-.l diffusivity 

B* B /B r 

V v' ref 

B„ Vertical conductivity ,y?C B„ 

' P V 

Cp Specific heat at constant pressure 
Euler number 
f. Coriolis parameter 

g Acceleration due to gravity 

h Depth relative to rigid lid 

H Reference depth, or vertical length scale 

I Grid index in x -direction or direction 

J Grid index in y -direction or .? -direction 

K Grid index in z-direction or ^-direction 

Kg Surface heat transfer coefficient 
L Horizontal length scale 

P Pressure 

P3 Surface or lid pressure 

Pj. Turbulent Prandtl number, 
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Pe I’fcU't nunilji.>r 

H. HuBiibv tminber 

Ue Heynoida number uarbuleuLj 

Hiclurdson number 
T Temperature 

Kevuretice lemperaiure 
Tjj Equilibrium cemperaturu 

T^ Surt'ace tumperaturt* 

t Time 

“ret time 

u veLv)ciLy in x-direction 

V VC loci: uV in v-dirccciun 

w veLoeiuV in n-diri*ct.ion 

X horir.unuai coordinite 

y horinontai coordinate 

c vertical coordinate 

Greek .^bmbol s 

t Horizontal coordinate in stretched system, ~ 

Horizontal coordinate in stretched system, » 

t Vertical coordi’aate in stretched system 

K/L 

Transformed vertical velocity 
Density 

r Surface shear stress in x-diroction 

vz 

f Surface shear stress in y-direction 
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Sup»?rsertp ti» 


( ) Dimensional 'maniiiy 

( ' ) Dimensional mean quantity 

( ' ) Dimensional fluctuating quantity 

( ) Dimensional quantity 

( ^g|leference quantity 


3 
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A]_ First terra in n^(t), which is deiined below 
A2 Coefficient of second term in 11 r^) 

3 j^ Horizontal eddy thermal diffusivity 

T1 

V Vertical eddy thermal r-iffusivity 

Cq Phase velocity or celerity of surface gravity waves 
Specific heat at constant pressure 
f Coriolis parameter 

g Accei.eration due to gravity 

h Depth relative to the mean water level 

H Depth contour relative to the free- surf ace, h + n 

I Grid index in x-direction or a direction 

J Crxd index in y-direction or 3 direction 

K Grid index in z-direction or a -direction 

Ky Horizontal kinematic eddy viscosity 

Vertical kinematic eddy viscosity 
K_ Surface heat transfer coefficient 

L Width O', bay at ocean-bay interface 

P Pressure 

P Surface pressure (atmospheric) 

s 

T Temperature 

T^ir Air temperature 

"^amb Water ambient temperature 

Tg, Equilibrium temperature 

Tg Surface temperature 

t Time 
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time lag in n^Ct), which is defined below 
u Velocity in x-direction (dimensional) 

V Velocity in y-direction (dimensional) 

Amplitude of inlet tidal velocity v^(t),VQ» 2A2 Cq 2 lll 

h ^ 

^o (t) Inlet tidal velocity ■ v^ cos aj(t+'J') 
w Velocity in z-direction (dimensional) 

X Horizontal coordinate 

y Horizontal coordinate 

z Vertical position relative to the mean xvater level 

Z Vertical position relative to the free surface, z -f- n 

Horizontal Stretching Parameners 
X Horizontal otretching coordinate in x-direction 

Y Horizontal stretching coordinate in y-direction 

da 

X" d^Y 

da^ 

y’ II 

de 

Y” d^Y 

a The distance at which minimum step size is desired in 

x-direction (see transformation relation below) 
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a,b , c-j_ , C2 . C3 c^ d and e are related as 
a = a + C]_ Sinn {c2 (X-d) } 

8 = b + c^ Sinh (c4(Y-e) } 


Greek Symbols 


a Horizoncai coordinate in stretched system, x 

3 Horizontal coordinate in stretched system, y 

•j Vertical coordinate in stretched system, Z/H 

U Transformed vertical velocity 

9 Density 

n Free-Surface elevation above mean water level 
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X Wave length of progressive wave at inlet 
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72 
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IWTRnDUCTION 


Background 

The management: of waste heat from power plants is a 
dominant consideration in making power production compatible 
with ecological concerns. For every unit of energy converted 
to electricity approximately two units are discharged into the 
environment. The ultimate heat sink is space. However, the 
intermediate heat transfer links, namely, the hydrosphere and 
the atmosphere may under, go changes harmful to life supporting 
ecosystems . 

Some quantitative estimates of efficiency, operating 
temperature and waste heat have been made by Harleman and 
Stolzenbach (1972) and are presented in Table I. Typical 
condenser water flow rate is about 1500 ft /sec (675,000 gal/min) 
or 3.4 X 10^16/hr. This results in about 12°F increase in 
cooling water temperature for fossil fuel and 20°F for nuclear 
plants. An idea of the magnitude of these discharges can be 
formed by observing the estimates given by Krenkel and Parker 
(1970) . According to their estimate the cooling water flow in 
the United States (based on a 15'^F rise in temperature) is 
approximately 40 trillion gallons per year which is approx- 
imately 10% of the total yearly flow of waters in the rivers 
and streams in the U.S. The problems are real. 

While the effects of thermal pollution have not been sys- 
tematically quantified, it is accepted that there are effects 
of significant nature in the biology and chemistry of the eco- 
system disturbed. Thermal discharges majr result in anomalous 
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•stratification of lakes, lowering of capacity to hold oxygen, 
increased reaction and increased metabolism. The lethal effects 
of thermal pollution are sometimes obvious, whereas the sub- 
lethal effects on food chains and waste assimilation capacities 
are not easy to foresee unless careful, fluid mechanical, 
chemical and biological interactive studies are conducted in 
an integrated fashion. 

In Che past, waste heat was primarily discharged directly 
by open cycle systems to aquatic ecosystems, eg. lakes, rivers, 
cooling ponds, etc. Recently, the shortage of land, particul- 
arly in Europe has resulted in closed systems, eg. cooling cowers 
chat discharge waste heat directly to the atmosphere. The 
incremental change in ecological impact implied by going from 
open to closed systems is still in the realm of investigation. 
Comparative statements are difficult to make especially when 
considering such non-ecological factors as economics . The pre- 
sent effort is directed solely towards hydrothermal analysis and 
predictions for open cycle cooling systems. 

Need for Models 

Accurate understanding of Che behavior of thermal discharges 
is important for the following reasons . 

I. To analyse the receiving body of water such that recircula- 
tion between intake and discharge from the condenser and con- 
sequent decrease in cooliig efficiency can be minimited. 

II. To assess Che thermal impact on the aquatic life forms 
existing in the receiving ecosystem, 

III . To provide a priori information about the nature and extent 

of thermal impact for site selection. 
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It is therefore, apparent that not only environmental but 
design Interests also are at stake. 

The above mentioned objectives can only be met by having 
large data sets over the entire discharge flowfield. .'jeasure- 
ments for temperature and velocity made over the affected domain 
could be used to develop maps for velocity and temperature 
distribution. However, there are some major drawbacks to this 
procedure . 

1. Unless the flowfield is adequately covered with fixed 
measuring installations which record temperature and velocity 
continuously, synoptic data is near impossible to obtain. In- 
situ measurements have been conventionally obtained from moving 
boats with towed measuring devices; the data consequently is 
non-synoptic giving distorted plume shapes except in rare sit- 
uations where a steady state plume exists. 

2. The information obtained is usually site and time specific 
and is quite difficult to use either for diagnostic or predictive 
purposes even at the same site under different meteorological 
and plant conditions or at other sites. Thus, the data obtained 
merely serves as a monitoring tool providing little information 
regarding the behavior of heated discharges that would be useful 
in analysis of existing discharges or design of outfall location 
for future power plants. 

3. Regulatory agencies often require data for worst meteorology 
situations which can only be obtained by simulation. 

4. For comparative studies of prospective sites and discharge 
geometry prior information regarding plume behavior is essential. 
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In-«ltu mcasurtnutiiLs are of no help In such sice selection 
decisions. 

While in>sicu measurements can serve for diagnostic and 
monitoring purposes under limited circumstances for .iteeting 
objective I and II they are not relevant for objective III. 

Models for simulating behavior of thermal discharges are there- 
fore imperative. 

Basic Considerations in Model Development 

In order to establish the rationale of model formulation 
the physical mechanisms underlying the heat dispersion from a 
heated discharge need to be outlined. Thermal discharges can 
either be in the form of surface jets or i the form of sub- 
merged jets. The surface jet is usually from a canal discharge 
whereas submerged jets are from either .single-port or multiport 
diffusers which release the thermal discharge at some finite depth 
below the air-water interface. Both for submerged jets and sur- 
face jets the following mechanisms govern the heat dispersal. 

1. Entrainment of ambient fluid into the thermal discharge. 

2. Buoyant spreading of discharged heated effluent. 

3. Diffusion by ambient turbulence. 

4. Interaction with ambient currents. 

5. Heat loss to the atmosphere through the air-water inter- 
face. 

The first four mechanisms redistribute heat and momenttim in the 
domain. The last mechanism transfers heat to the atmosphere. 
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Two factors which affect plume behavior at'e 

1. Dlscharfc,e geometry and location with respect to 

ambient stratification. 

2. Interaction of discharge with bottom topography. 

The mef mnlsms mentioned play roles of varying importance 
as the heated effluent travels away from the mouth of the di.s- 
charge. It has been customary therefore, to divide the flow 
into the following regions; 

1. :?ear>Fteld 

In this region the Initial properties of the discharges 
are imoortant. The flow field Is dominated by the jet like 
structure of the discharge. Discharge geometry, bottom topo- 
graphy, initial velocity, temperature, etc are dominating 
variables, Thus, non-dimensional qualities such at jet 
Reynolds No., densimetric Froude A). , aspect ratio, bottom slope, 
anu jeu uepth to domain depth ratio are important. Redis- 
tribution of heat from the discharge to the ambient is pri- 
marily by entrainment of ambient fluid. Non-significant amounts 
of heat is transferred to the atmosphere. 

2 . Far-Field 

Here the ambient conditions are dominant in heat dxsper.sal. 
Ambient turbulent diffusion and surface heat loss are signif- 
icant. 

The boundary between near and far field is quite quali- 
tative with no easily definable senaration 

Owing to the relative importance o; different heat and 
mass transfer processes in the two regimes and consequent 
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different approxlTmtlon* it has been customary to develop 
different models for each of the domains. Models that con- 
sider the complete domain are termed complin e-fleld models. 

It is Important to note that almost all models to date do not 
simulate the ambient conditions but input them as boundary 
conditions obtained by measurement. 

The approach in the models developed in the present effort 
has been to develop complete field models, which not only 
simulate the thermal anomaly region, but also the total 
ambient condition. Therefore, the mcdeis are comprehensive 
and the following definitions of near and far fields are 
stipulated. 

Near-Field - The complete region where effects of thermal 
discharges show measureabie and distinct thermal perturba- 
tion to ambient conditions. Thus, this definition includes 
near, and far fields of traditional definition. 

“• Far- Fie Id - This is the total domain •..r ose thermcHlynamic -md 
hydrodynamic characteristics affects the discharge, but is large 
enough not to be significantly affected by the discharge. 

Thus the far-field solution affects the near field boun- 

» 

dary conditions, whereas the far-fieid is significantly un- 
affected by the near field characteristics. Throughout the 
three volume report these definitions of near and far field are 
used. 

The first step for model development is to mathematically 
represent the conservation of mass, momentum and heat in terms 
of a set of equations. A relationship between the temperature 
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and velocity field in the form of an equation of state completes 
the set of equations. The second step involves assumptions and 
approximations for different flow regions in order to make the 
equation mathematically tractable. The final step is to 
develop a solution procedure which obtains solutions with 
appropria a boundary conditions . 

Review of Other Modeling Efforts 

An excellent review and evaluation of 40 sutfacea plume 
models have been presented by Dunn et al (1975) . Salient 
featuies of past modeling efforts will be highlighted in this 
section. 

Models can be classified in the following groups: 

1. Phenomenological Models 

These models are basically empirical correlations of num- 
erous data bases. Measured plume characteristics such as 
centerline temperature decay, jet width, isotherm areas are 
nrelated with jet and domain variables such as initial 
^ensimetric Froude No., Reynolds No., bottom slope, outfall 
geometric parameters, etc. These models are relatively easy 
to use. However, information regarding detailed distributions 
are not available from such models. Another major disadvantage 
is that these models are only applicable to specific physical 
situations for which correlations were obtained. The models 
developed by Carter and Regier (1974) and Pritchard (1973) 
are representative examples of such models. While these models 
give gross parametric descriptions of plume behavior they are 
not useful in analysis of recirculation, interaction with 
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ambient currents .winds , tides, etc., or in understanding the 
time dependent dispersion of heat in the receiving water body 
under varying boundary conditions. 

Integral Models 

In these models only regions arbritarily defined as the 
pliime are considered. Forms of velocity and temperature pro- 
files normal to the axis of the jet are assumed to provide 
cLoaui.e for integrated conservation equations. Numerous 
models of this type are in existence with varying degrees of 
simulation success depending on discharge and domain geometry. 
There are some basic deficience.s in these NoHels. 

1. Domain boundaries art not considered. 

a) Domains are considered to ha>re sufficient depths to 
eliminate bottom boundary flow effects. 

b) Effects of lateral boundaries are completely neglected. 

2. Ambient stratification is ignored thus resulting in serious 
errors in incorporating bouyancy effects. This is a major 
limiting feature. While integral models may sometimes be ade- 
quate for non-bouyant jets, their applicability for bouyant 
jets is almost universally questionable. 

3. Changing ambient currents can not be incorporated. 

4. Wind effects are ignored. 

5. These are steady state models, and verification is near 
impossible in field situations owing to near impossibility of 
encountering steady state plumes. 

6. Entrainment coefficients are a function of numerous jet and 
.'’™hient parameters making generally acceptable coefficients 
difficult to compute. 
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This class of models has been widely used since they have 
some predictive capability and are computationally economical 
compared to numerical models. Table II shows a list of re- 
presentative integral models. Sharazi and Davis (1974) have 
developed workbooks using Prych (1973) integral type model. 

They present numerous nomograms facilitating use of this model. 
Numerical Models 

The state at a point in a flow field la described by the 
solution of a system of equations that describe the local 
conservation laws for total mass, species mass, momentum and 
energy. For thermal pollution modeling the conservation laws 
for mass, momentum and energy are relevant. The constitutive 
equations complete the set. Since most environmental flows are 
turbulent a closure condition for the turbulence model is re- 
quired. This system of equations together with appropriate 
boundary conditions constitute a mathematical model. The 
equations are coupled, non-linear, second order, three dimen- 
sional partial differential equations . Analytical solutions are 
not possible for most practical situations. Various assumptions 
regarding dimensions and relative importance of physical mechanisms 
are required to make the equations mathematically tractable. 

The differential numerical approach attempts to find solutions 
to this system using numerical techniques. With the advent of 
high speed computers and appropriate numerical techniques this 
approach in thermal pollution modeling has become increasingly 
popular. The main advantage of this procedure compared to the 
others described in previous sections is in the promise of ade- 
quate simulation of all important physical mechanisms without 
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Che need for damaging assumptions regarding the nature of the 
flow. However, assumptions and approximations can still be 
made for specific situations. More importantly the three-di- 
mensional nature of bouyant plumes can be accurately simulated 
by this approach. These models also have the capability of 
simulating time dependent behavior with time varying boundary 
conditions . 

A number of models of this class are in existence. Table 
III shows a representative group of thef.e models. It can be 
seen that though there are a number of three-dimensional models 
available all have some limiting assumptions. The purpose of 
the present modeling effort was to develop a model package that 
could be applied to a large variety of discharge, domain and 
meteorological conditions with relative ease. The existing models 
mainly suffer for inadequate ability to include bottom topography • 
effects, and ambient meterological conditions. They also lack 
adequc verification. 

Rationale for Present Models 

A report by Lee and Sengupta (1976) presents the details of 
model development. The present section presents the summary 
of this effort. 

The Thermal Pollution research team at the University of 
Miami, under the sponsorship of NASA-KSC, has been for the last 
few years developing a package of mathematical models which can 
have general application to problems of power plant heated 
discharge to the aquatic ecosystem. The effort is closely in- 
tegrated with simultaneous remote sensing and ground- truth 
data acquisition support. The concept being the development of 
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adequately calibrated and verified models for direct application 
by the user community. The user community being the utilities 
and the regulatory agencies like the Environmental Protection 
Agency and the Nuclear Regulatory Commission. 

Critical evaluation of mathematical models in use for 
thermal pollution analysis, indicated that though some models 
may perform well under certain conditions , a generalized three 
dimensional model which accounts for wind, current, tide, bot- 
tom topography and diverse meteorological conditions was not 
in existence. The NASA-KSC-University of Miami project has 
specifically proceeded to develop a model package which satis- 
fies these requirements with emphasis on remote sensing data 
input and verification during the model development. The com- 
plete effort in flow chart form is shown in Fig. 1-1 to indicate 
the relationship between data acquisition and model development. 
The Model Package 

Critical evaluation of mathematical models used for thermal 
pollution analysis has been made by Dunn et al (1975) . They 
compared the performance of various models in predicting a stan- 
dard data base. A general conclusion that can be made from their 
analysis is that though some models may perform well under cer- 
tain conditions a generalized model which accounts for wind, 
current, tide, bottom topography and diverse meteorological 
conditions is yet to be developed. The models in existence can 
be classified into three categories: integral models which 
make similarity assumptions, phenomenological models which rely 
heavily on data and numerical models which solve boundary value 
formulations directly using numerical tehcniques. Finite 
difference methods are widely used though some finite-element 
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algorithms have been tried. The integral models of Motz and 
Benedict (1972), Prych (1972), Shirazi and Davis (1974), and 
Stolzenbach and Harleman (1972), are the ones that have been 
applied to a number of situations with mixed success . The two 
phenomenological models that have been widely used are by 
Pritchard (1971, 1974). The numerical models can be classified 
in terms of spatial dimensions used and physical assumptions 
made. There are a number of two dimensional models that have 
been developed by Trent (1974, 1975). The formulations even 
for the three dimensional methods vary widely from primitive 
variable to velocity-vorticity and velocity-corrector potential 
methods. One of the first three dimensional models was by Waldrop 
and Farmer (1973, 1974, a,b). This model was essentially a 
free surface formulation. One of the first three-dimensional 
models which adequately accounted for bottom topography and 
comprehensive meteorological conditions was a rigid-lid model 
developed by Sengupta and Lick (1974a, 1976b) . They used a 
vertical stretching to convert a variable depth basin to constant 
depth. Irregular shorelines could be easily included without 
modification of the program. Modified versions of this model has 
been used at a nximber of sites with satisfactory results, e.g., 
Sengupta and Lee (1976a) , as part of the University of Miami 
model development effort. 

Most mathematical models for environmental flows require 
detailed verification. The nature of the governing equations and 
the state of the art in solution techniques demand restrictive 
assumptions and approximations. Often boundary conditions and 
initial conditions are not adequate. Thus a careful calibration 
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procedure is an integral part of a model development effort. 

It is not sufficient to verify the model after it has been 
developed. Simultaneous calibration as the model is developed 
leads to modifications which check whether assumptions are 
valid and may even help to simplify models. A detailed dis- 
cussion of the need and methodology for calibration and ver- 
ification of ntimerical models is presented by Sengupta et al 
(1975) . The present study incoporates a program of airborne 
radiometer data and in-situ measurements to enhance the model 
development effort. Figure 2-1 shows the interrelationship of 
the data gathering and model development efforts . The details 
of the mathematical package and formulation are presented in a 
number of reports by Veziroglf et al (1974, 1975), and a summary 
presentation is given by Sengupta and Lee (1976a) . A brief 
discussion will be given here. 

The primary motivation behind the effort was to develop 
a series of models which make minimal site restrictive assumptions 
enabling application to diverse basin and discharge configur- 
ations. Two separate formulations were made one with the rigid- 
lid approximation and the other with the free-surface included. 

The rigid-lid formulation was essentially an extension of 
the effort by Sengupta (1976b), to facilitate application to 
thermal pollution studies. The free-surface formulation is 
similar to that of Freeman et al (1972) , used in a study of Lake 
Huron. The models are further modified to have specific app- 
lication to near field and far field. The near field being that 
region affected by the pltune and the far-field being the encom- 
passing domain. The far field affects the near field. The near 
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field has minor effects on the far field. The near field models 
are especially capable of modeling open-boundary conditions. 

Thus there are four separate versions of the program, near and 
far field versions of rigid-lid and free-surface models. All 
the models include a vertical normalization with respect to 
local height to convert variable depth domains to constant 
depth. The programs have a horizontal grid-point marking 
system which allows application to different shore line geometry 
without any modification to the program. One version of the 
free-surface model has a hyperbolic sine stretching similar to 
Waldrop and Farmer's (1974a), tangent str'»tching to allow finer 
resolution at discharge points. Figure 1-2 shows the component 
programs of the mathematical package and present application sites. 

The governing equations for the rigid-lid model are the 
Incompressible Navier-Stokes Equation, conservation of mass, 
energy and an equation of state, A predictive equation for lid 
pressure is derived from the vertically integrated horizontal 
momentum equations, The hydrostatic, Boussinesq and rigid-lid 
approximations are made. The turbulent closure condition is made 
by using eddy transport coefficients . The boundary conditions 
at solid boundaries are no-slip, no normal velocity and adiabatic 
conditions, At the air-water interface wind stress and heat 
transfer coefficients are specified (a conduction formulation) . 

At open boundaries conditions are specified for temperature and 
velocities where available. Otherwise normal derivatives are 
equal to zero. Complete conditions in space and time are 
specified at discharge locations. Explicit numerical schemes 
with Dufort-Frankel differencing of the diffusion terms are used. 
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The pressure equation is solved by the SOR technique. 

The free surface model also uses the same equations except 
the predictive equation for pressure is the hydrostatic equation. 

The surface pressure is atmospheric. One extra equation for sur- 
face height obtained by vertically integrating the continuity 
equation is used. The other approximations are the same as that 
for the rigid lid formulation. The boundary conditions are also 
the same except that at open boundaries conditions for surface 
height are required especially for tidal situations. Lateral 
walls have slip conditions. 

The mathematical model package therefore consists of: 

1. Rigid-Lid Model 

a) Far field version 

b) Near field version 

2. Free surface model 

a) Far field version 

b) Near field version with horizontal stretching 
The features of the models can be detailed as follox/s : 

Features of rigid- lid model 

Thr e e - dimen s iona 1 
Non-linear 
Baroclinic 
Time -dependent 
Irregular topography 

Driving forces: wind, heat and mass flux 

Predicts three-dimensional fields for velocity and temperature 
Surface pressure is defined as the pressure on a rigid-lid 
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Features of the free-surface model 


Three-dimensional 
Non-linear 
Baroclinic 
Time-dependent 
Irregular topography 

Driving forces: wind, tide, heat, and mass flux 

Predicts three-dimensional fields for velocity and 
temperature . 

Predicts surface height. 

Application Sites 

The application sites were chosen to represent as diverse 
a group of topographical situations as possible. The sites 
were Biscayne Bay in South Florida, Hutchinson Island in mid- 
Florida along the Atlantic coast and Lake Be lews in North 
Carolina. Biscayne Bay shown in Fig. 1-3 is a shallow estuarine 
basin with tidal exchange with the Atlantic Ocean through a 
safety valve region and a number of creeks . There are two 
power plants on the Bay, operated by the Florida Pov;er and 
Light Company. The Cutler Ridge Plant shown in Fig. 1-4 
has a canal surface discharge into the Bay. The Turkey Point 
Plant has a closed canal cooling system. .Florida Power and Light 
Company has a newly built plant at Hutchinson Island. This 
is an open ocean, coastal discharge approximately 1200 feet 
offshore. The discharge mouth is Y-shaped pipe with one leg 
60° to the other. This is a submerged discharge. Fig. 1-5 
shows the site. Lake Belews of Duke Povjer Company in North 
Carolina has a mixing pond for the heated discharge. Connected 
by a canal to a larger cooling lake where plant intake is 
situated. This is a t ’pical man-made lake for the Southeastern 
United States. Fig. 1-6 shows the lake configuration. The 
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nixing pond ia themaily well mixed whereas the main lake dis- 
plays a seasonal thermocline. This is a surface discharge 
situation. 

Detailed results for each application site are presented 
in reports by Lee et al (1974, 1975, 1976) and in dissertations 
by Venkata (1977) , Mathavan (1977) and Tsai (1977) Applications 
to each site was verified with a number of data bases collected 
during field experiments using boats and fixed stations for 
ground truth and in- situ measurements and airborne radio- 
meter for remote measurements of surface temperature. Cali- 
bration and verification was conducted by specifying detailed 
initial conditions and comparing simulated results for given 
time intervals with subsequent data bases. Thus, the capability 
of the models to include varying boundary conditions in sim- 
ulating time dependent behavior was severely tested. 

Conclusions 

The salient conclusions of thu project, to date, are as 
follows : 

1. The importance and need for airborne remote sensing data 
has been demonstrated for thermal pollution studies. 

2. The unique role of synoptic data bases obtained by remote 
sensing, has been clearly established as imperative for 
complex model development efforts, 

3. A mathematical model package has been developed with 
adequate inclusion of complex transport processes to serve 
as a predict! ’"e tool for thermal discharge studies and 
site selection. 

4. The component programs of the package have been applied to 
diverse discharge and topographical conditions . The models 
have performed satisfactorily for different meteorological 
inputs , 
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5. It can be concluded that a reasonably general model package 
has been developed for applii-acion by the user cc..jriunity . 
Minimal programming effort is required by the prospective 
users of the model package. 


Recommendations 

I-Jhile the models have been verified for a number of sites, 
in order to inspire greater user confidence it is Imperative 
that applications to at least two other sites be made. One 
site should be such chat the rigid-lid model is appropriate. 

The other site should have the features which test the cap- 
abilities of the free-surface model. Careful verification of 
velocity prediction should also be emp'^isiaed. 
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TABIS HOt II LIST OF THE INTEGRAL LEVEIiOPED 
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MA TimMATlCAl. MODEL PACKAGE 


RIGID LID MODEL 


FREE- SURFACE MODEL 



Par Field Version 

i. 

Far Field Version 

■ 

Biscayne Bay 


Biscayne Bay 

ii. 

Near Field Version 
cutler Ridge Plume 

f 

ii. 

.Near Field Version 
With Horizontal 
Streching 
Hutchinson Island 

iii. 

Verification Site 
Belews Lake 




ALL MODELS INCORPORATED A NORMALIZATION WITH RESPECT TO LOCAL DEPTH 


Figure 1-2 
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Fig. 1-3 Map of southeastern Dade County showing the area 
of investigation and location of gaging stations. 35 
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* Th«.Rlqld-Lid Model 

2,1 Brl«f Dggcriptlon of Pa«t Exparlence with Rigid-Lid 
Fonmilationt T Relevant Advantages and Diaadvantaqes 
Ona of tha first bhraa-dlmansional modals which adaguataly 
accountad for bottom topography and comprehansiva mateorological 
conditions was a rigid-lid modal davalopad by Sangupta and Lick 
(1974, 1976) . Thay usad a vertical stretching to convert a 
variable depth basin to constant depth, thereby permitting a 
constant vertical grid size to be applied everywhere in the 
domain, irregular shorelines could be easily included without 
modification of the computer program. Modified versions of the 
rigid-lid model have been used, by the Thermal pollution Re- 
search team at the University of Miami, at a number of sites 
with satisfactory results; Lee and Sengupta (1976), 

The major advantage associated with the rigid-lid model is 
the elimination of surface gravity waves with a consequent larger 
integration time steps. For sites where gravity waves do not 
determine the maximum allowable time step (for example, in the 
case for which vertical diffusion determines the maximum allow- 
able time step. This could be the case for a shallow water 
basin), no time step advantage is gained by the rigid-lid model. 
The major disadvantage associated with the rigid-lid model 
is its inability to predict surface heights. 

Thus, for example, real tidal conditions cannot be accounted 
for, since the surface is not permitted to move, 

2.2. Assumptions and Approximations 

The system of governing equations (see next section, 2.3) 
for the fluid flow invoke several simplifying assumptions and 
approximations in the interest of saving computational time 
without losing significant accuracy. The following assumptions 
and approximations have been employed: 
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2.2.1 The Bousslnesq Approxlmaelon 

The effect of density variation on the inertial and 
diffusion terms in the governing conservation equations is 
neglected. Density variation is retained in the bouyancy terms 
in the equations of motion. The effect of bouyancy is thereby 
accounted for by allowing density variations which produce 
horizontal pressure gradients which influence the fluid motion 
through the horizontal momentum equations. 

2.2.2 The Hydrostatic Approximation 

The hydrostatic approximation involves neglecting the 
vertical convection and diffusion terms in the vertical momentum 
equation. This approximation implies that the vertical fluid 
acceleration, is negligible. 

2.2.3 Constant Eddy Transport Coefficients 

Turbulence modeling is very complex and has an extensive 

body of literature of its own. Turbulent closure has been 
obtained in this model by using constant eddy transport co- 
efficients, except for the case (Lake Belews site) for thermal 
stratification. For this case a Richardson number dependent 
variable vertical eddy transport coefficient was used. Due 
to the horizontal scale length, L, being much larger than the 
vertical scale length H, the horizontal eddy transport co- 
efficient is orders of magnitude larger than the vertical eddy 
transport coefficient. 

2.2.4 Variation of Surface Wind Stresses 

The variation of the wind produced surface shear stresses 
with respect to x and y, and , are considered 

negligible for the horizontal length scales of the water bodies 
studied. However, if the physical dimensions of the water body 
are so large as to require including variation of the wind 
stresses with respect to x and y, then the computer programs 
can be quite easily modified by replacing and with 
matrices and Ty.^(I,J). Where the indices I and J 

refer to the location of a grid point with respect to the x,y .plane 
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2.2.5 The RtRld~Lid Approximation 

This approximation effectively eliminates surface gravity 
waves by imposing a zero vertical velocity at the surface. 

This means that the surface allows slip conditions without any 
normal velocity. In other words, a rigid, ^^frictionless lid has 
been placed at the undisturbed free surface of the water body. 

The surface pressure is no longer atmospheric, but represents 
a "lid-pressure", which under some special steady state conditions 
can be related to the free surface elevation that would occur 
if no lid were present. 

2.3 Governing Equations 

The set of equations governing the behavior of the fluid 
flow are those expressing the conservation of mass , momentxam 
and energy in turbulent flow, and an equation of state. 

2.3.1 Cartreslan Coordinate Representation (x,y,z) 

The Cartesian coordinate system is used with the z- 
coordinate in the downward vertical direction as shown in Fig, 

2.1, i.e. a so-called "left-handed" coordinate system. In 
order to keep the generalized nature of the model, all the 
significant terms in the respective conservation equations are 
retained. Included are the effects of bouyancy, inertia, 

Coriolis, density and turbulent mixing. Wind shear and heat 
flux at the surface are also considered. 

The following system of non-linear partial differential 
equations, written in cartesian coordinates , describes the 
three-dimensional: unsteady fluid flow where the variables are 

in non-dimensional form. 

Continuity Equation 

3u + 3v + 3w “0 (2.1) 

3x 3y 3 s 
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Horizontal ‘.'Momentum Equations 

+ u3 U + VUl + W )U - I 

3t ix 3y jzRij 
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(2.3) 


Hydrostatic Pressure Equation 
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E^(l+ p) 
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Energy Equation 
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(2.5) 


Equation of State 


P * P ( ^) (2.6) 

Where the set of non-dimensional quantities are defined as 
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where quantities widh the subscript 'ref are reference quantities 
for the respective variables. The tilda denotes dimensional 
quantities. The quantities with an asterisk superscript are 
non-dimensional quantities. 

L is the horizontal length scale. 

H is the vertical length scale. 

The other symbols in equations (2.1) - (2.6) are defined in the 
list of symbols for the rigid-lid model. 

2.3.2 Vertical Stretched Coordinate Representation C^^ ^ ^ ■:) 

The programming difficulties for a three-dimensional 
basin suggest a stretching of the vertical coordinate of the 
form 

- z (2.70) 

(x.y) 

This coordinate transformation converts the basin to a constant 
depth one, so that a constant vertical grid size, aV , can be 
used throughout the domain. The horizontal coordinates (x,y) 
are transformed by letting 

^ ^ ( 2 . 8 ) 

2 - y 

Fig. 2.2 shows the coordinate system for 100 grid points . (sample) 

Once again although the symbols used in this section are defined 
in the list of symbols for the rigid-lid model, it is worthwhile 
to explain here the meaning of the subscripts and wavy lines 
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(cildas) used In the following vertically stretched equations. 
Quantities with subscript 'ref are reference quantities, H 
and L are vertical and horizontal length scales, respectively. 
The variables with wavy lines on top are dimensional quantities. 
The quantities with an asterisk superscript are non-dimensional. 

By converting the basin to constant depth, the same number 
of grid points, and of constant grid size, can be used in the 
vertical direction in shallow as well as in deep regions . The 
details of transforming the equations in the (x,y,z) coordinate 
system into the coordinate system is given in Ser.gupta 

and Lick (1974) 

The non-dimensional governing equations in the ^ 

coordinate system expressed as follows: 

Continuity Equation 

jihiil + h|| -0 (2.9) 

Horizontal Momentum Equations 
(hu) ^ 3(huu ) + iLhwd +h^j^ - h 

. 3 Ps - h + 1 1_ (hDu) + 1 :• (h3u) 

3 1 ( KT t 
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( 2 . 10 ) 
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where 

3P - B + 3Ps 
^ 3 X 


3P - + 21i 

3y ' 


and 

43C 0 


B„ » E,, 3 ■ odz 

^ ^ Ty o 


and 

« Surface Pressure 
s 


Hydrostatic Equation 


DP - E, (H-p)h 
3? ^ 


( 2 . 12 ) 


Energy Equation 


9(hT) a (huT) . \(hvT) + h a(;7r) 
3 1 3 u 3 ^ 9 i 


1 

P7 


e 


1 (h +L. 1_ rhil'i + 1 1 3 fv,* 

3a 3a^ ^ 33 ^*^33^ PTi^ E ^ 


(2.13) 

Equation o£ State 

p « p (^) is given for fresh and salt water as follows: 

Salt Water: p (^) » 1.029431 -.000020^ -.0000048^ ^ (2.14) 

(for salinity of 38 parts per thousand^ 

Fresh Water :p ('i?) » 0.000428 -.000019^ -.0000046't^ (2.15) 


where again the wavy line denotes dimensional quantities. Note, 
^ is in degrees Celsius . 

Actual Vertical Velocity in ((K,'y,'z) coordinate system 



(note, by virtue of the rigid-lid approximation, w (z»=0)=0) 
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Due to the rigid-lid approximation the surface pressure, 

P„, is not atmospheric as in the case for a free surface model, 
s 

To obtain a predictive equation for surface pressure, P , the 

s 

horizontal momentum equations (2.10) and (2.11) are integrated 
from z»o to z«h, where h is the non-dimensional depth h/H. 

The integrated equations are then differentiated once with 
respect to a and 3 and then summed. This derivation yields the 
Poisson equation for surface pressure, P^ . Sengupta and 
Lick (1974). 

Surface Pressure Equation 
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(2.17) 


The last term is the Hirt and Harlow (1964) correction 
term which accounts for non-zero vertical velocities at the 
rigid-lid. The variables (b^, B y , and \l. A ^2 etc) are given 
below: 

"^xi */ h 

0 

Av 9 “ ^ vd2f 

X.- f 
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3 (h3u) + 3_ (h3u) + Ip 1 3_ (A‘^ 3 u) 7 d 

3a 3a 3 3 33 H 5 y TT 
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2.4 Boundary Conditions (closed basins and open basins) 

The nature of the system of governing equations requires 
initial and boundary conditions to be specified. The boundary 
conditions for both near-field (open basin) and far-field (closed 
basin) versions of the rigid-lid model are presented in this 
section. The initial conditions will be presented in the next 
section. 
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2.4.1 Near-Field (open bastu) 

The rigid-lid model, near-field has been applied to an 
open boundary domain with thermal discharge along the lateral 
solid boundary, Venkata (1977). This model has also been applied 
by Mathavan (1977) to a mixing pond which has an opening at 
one boundary and a thermal discharge at another boundary. 

The set of boundary conditions for the domain having open 
boundaries will be given in this section, since these open 
boundary conditions are most difficult to specify. 

At solid boundaries no-slip and zero normal velocity conditions 
are specified. 

At the air-water interface 


All solid boundaries are considered adiabatic. 

1^ and |~ are set proportional to 

the wind shear stress components in the x and y directions, 

respectively. The rigid-lid approximation sets the vertical 

velocity w to be zero at the air-water interface. Also at the 
3T 

surface is set proportional to the heat flxix from the 
surface. The heat flux from the surface is in turn proportional 


to (T^-T^), 


whore T^ is the so-called equilibrium temperature, 


The equilibrium temperature is the surface temperature at which 
the heat entering the water body at the air-water interface is 
equal to heat leaving the water surface. At the open boundaries 
the first order derivatives of temperature and velocity are 
set equal to zero. Thus, the boundary conditions are in summary: 


At the surface, 

J2 =* 0 

( 1 ) 


At lateral solid boundaries 
on x-boundaries : 


3u 

- -( _ JlH i 

“ 0 

3<r 

“refP ^ 
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3v 

= - ( hH )Tzy 

V =■ 0 


Uref^\ 

3T 

3/ 

= - ( hHKs )(T -TJ 
^ e s 

3T = 
3x 


3tjt 


- X 
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3h 
3 a 


oT 

3^ 


( 1 ) 


Note: H is vertical scale length for rigid-lid model. 
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At the bottom. 
i2 - 0 
u 0 

V » 0 

3T - 0 

At lateral open boundaries 
on x-bovindarles 

- 0 

( 2 ) 0 

3a 

V - 0 

DT .3T_Y^3T»0 
3 X 3 a n 3 a 3 

At Discharge 

Velocity ? 

Density ( Specified 

Temperature/ 

2.4.2 Far -Field (closed basin) 

The rigid-lid model, far-field has been applied to the 
Biscayne Bay, Sengupta (1975), and to Lake Belews, Mathavan 
(1977) . For both applications there is no direct thermal 
discharge modeled, but open boundaries are treated in much 
the same manner as was outlined for the near-field studies. 
The boundary conditions are in summary: 


At the surface, ^ 

( *0 

At 

lateral solid boundaries 

-Q » 0 


on 

x-boundaries 

3u = -( hH ) 

TZX 
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^ “re£»\ 
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3v ** - ( hH ) 
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= ^ - r 3h ")T = 0 


3x 
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8T - -( hHKs ) 

(T -T ) 



T? PCpB^ 

e s 




Note: If u is into domain then u=u far field value. If 
u points out of domain then u . j adjacent interior 

point. Similarly for vspecificatrony^ 
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o n y-boundaries 

ii • 0 

u - 0 
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3T-3T - ^ 3h 3T -0 
3y 3B h 3B at 
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2.5 Initial Conditions 

The initial conditions are specified by using the 
corrected IR data base for temperature; as will be 
illustrated in the sample problems in Volume II and zero 
velocity eT^erywhere in the domain (u-iva.Si-0 ), since it is 
quite difficult to obtain ground truth current measurements 
for the entire domain for the kind of grid size resolution 
that would be required. 

2. 6 Method of Solution 

2.6.1 General 

It is obvious that closed- form analytical solution of 
the system of governing aquations (2.9) - (2.17) is impossible 
to get. The set of equations consists of coupled, unsteady, 
three-dimensional, nonlinear partial differential equations. 
Therefore, the finite difference method is used to obtain 
ntmerical solutions. 

A three-dimensional grid system is established with 
respect to the ( a, 8,Y) coordinate system for the rigid- 
lid model. The governing equations are then solved over 
finite time steps which are carefully selected to obey nxom- 
erical stability criteria. This will be discussed in detail 
in a following subsection on stability criteria (2.6.5). 

^^^Note: At an inlet either u(t) or v(t) must specified, and at 
an outlet ^ specified. 

3a 96 
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In general, several methods are available for integ- 
rating over time the governing equations for incompressible 
fluid flow as discussed by Roache (1972) . The two most 
common techniques for integrating time-dependent partial 
differential equations are the implicit and explicit finite 
difference methods. The implicit method involves the sol- 
ution of a set of simultaneous equations which are obtained 
by writing the spatial derivatives in teriiis of the res- 
pective unknown quantities at the current time level, n+1, 
knowing the values of the remaining quantities of the set 
( u,v,JJ.P»P»T ) at the two previous time levels n and n .1, 

For the one-dimensional case, implicit methods are convenient 
because the set of simultaneous equations is tri- 
diagonal, Richtmyer and Morton (1967), and, hence, a direct 
matrix inversion method of solution is used. However, 
in the case of a three-dimensional model, the implicit 
method becomes too time constaming; since the simultaneous 
equations must be solved at each time step by an iterative 
technique. Thus, although the advantage of implicit methods 
is that they allow larger time steps, for the three dimen- 
sional case the iteration time for each time step more than 
offsets the inherently larger time step. Furthermore, 
alternating direction- implicit (ADI) methods may be used 
to obtain tri-diagonal matrices even for multidimensional 
equations, however, for irregular boundaries the ADI methods 
are impractical. 

Therefore, the explicit finite difference method is 
used for numerical solution of the rigid-lid model. The 
solution to a particular partial differential equation is 
propagated from point to point on the numerical grid system. 
The current time level value n + 1, of a particular system 
variable (u,v/l,P, p,T) is computed in general from known 
values of the corresponding system variables at the two 
previous time levels n and n-1. Thus, this is an explicit 
scheme. However, as will be seen later, the governing 
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equations for surface (or lid) pressure for the rigid-lid 
modal is elliptic; Sengupca and Lick (1974), and, therefore, 
values of for the entire domain are computed at each time 
step, iteratively. 

The mathematical model is an initial-value, boundary 
value problem and, hence, requires specification of both 
initial conditions and boundary conditions (see sections 
2.4 and 2.5) . 

2.6.2 Computational Grid System 

A schematic of the computational grid system for the 
rigid-lid model is illustrated in Pig. 2.3. A half-grid 
(dashed-line grid) is superimposed on a full grid (solid 
line grid) in the ag plane. The horizontal velocity 
components u and v are defined at the nodes of the full grid 
at (I,J,K); and P,P, T are defined at the nodes of the 
half-grid ( K). This arrangement allows better 
meshing of the system variables at all the nodes of this 
staggered mesh system. Constant grid spacing is used in the 
Y direction. Constant grid spacing may or may not be used 
in the a and g directions, respectively. 

During computation) values of the system variables 
u and V in the full grid system (I,J,K) are averaged, 
as follows , for computation of the system variables iJ.P.p 
and T in the half-grid system J-f*’.;, K) : 


(u,v) 


I+%, J+?2 


C(u,v) 


I.J.K 


+ (u,v) 


I+l, J+1,K 


+(u,v) 


I, J, -11,1: 


2.6.3 MAR and HRtl Numbering System 
The computational full grid system is divided into 
separate region:-? depending on the type of spatial finite 
difference used. That is, a tV 70 -dimensional matrix call 
MAR (I,J) is used in the model v;hich distinguishes between 
interior points, points on the boundary j and points out- 
side the domain of solution. Similarly, for the half-grid 
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system a two-dimensional matrix call MRH ( I + J + 
is used to distinguish between spatial finite differencing 
in the interior, on the boundary, and outside the domain of 
solution; Sengupta and Lick (1974). The MAR numbering system 
and the MRH numbering system will be clearly specified in 
the sample problems in Volume II. 

2.6.4 Finite Difference Schemes 

2. 6.^.1 Approximation of Spatial and Temporal Partial 
DerivativesjConservative Form 

The spatial derivatives are central differenced in the 
interior ; for example in the full grid system; 


3u ^u(I+l, J. K) -U(I-1. J, K) 
la JIa 

and, 

1 ud+l.J.K) +uq-l.J.K^ -2ua.J.K) 

and in the half-grid system: 

II 2; T(I+ 3/2. J.K) -T(I-%.J.K) 

3a TEa 


(2.19) 

( 2 . 20 ) 


( 2 . 21 ) 


and, 

^ i T(I+3/2,J,K) + ig-ij.J.K) - 2T(H-V.J.K) (2.22) 

^ (ia)=' 

At the boundaries . however, three-point single sided schemes 
are used by fitting a parabola through three points (the 
boundary point and the next twc coincident interior points) . 


Thus, for example, at the left a- boundary: 

‘X, 


3u 

- 4u(I+l,J,K) -3u(I,J.K) -u(I+2.J,K) 


3a 

2Aa 

(2. .23) 

and, 



3^u 

— S’ 
3a 

2; u(I,J,K) +u(I+2.J,K) -2u(I+l,J,K) 
(Aa) 

(2.24) 

and, 

at the right a -boundary: 


3u 

= 3u(I,J,K) +u(I-2,J,K) -4u(I-l.J.K) 


3a 

2Aa 

(2.25) 
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(2.26) 


and, 


— I 
da^ 


u(I,J,K) •Hag-2, J.K) >2 u(I-1,J.K) 


Similar txprassion to aquations (2.23) - (2.26) may ba 
obtainad for systam variablas in tha half -grid syatam by 
simply raplacing (I+1,J,K) with (1+3/2, J.IJ, ( I-l, J,K ) 
with (I-%,J,K ), (I, J, K. ) with (I+%,J,K ), (I+2.J, K ) 

with ( 1+5/2, J,K), and (1-2, J, K.) with (I-3/2,J,K). 

Nota, that tha spatial finita difference approximations 
(2.19) -(2.26) are on order of accuracy of ( ia )^. Crandall (1965). 

The temporal derivatives can be expressed in two forms, first 
n+1 n 

^ ^ u(I.J.K) - uCI.J.K) (2 ..27) 


3t 


1C 


for forward differencing in time) which is on the order of 


accuracy of and 

TH ' iic ( 2 „ 28 ) 


for central differencing in time; which is on the order of 

2 

accuracy of( At) . 

The finite differences in both space and time in the 
model are expressed in the full conservation forms following 
Arakawa (1966), for example: 

3 (Hu) - ^^^h+l,J,K 1-1, J.K C2 s29) 

9 a ila 


This is done to avoid possible "leaking" of mass, momentum, 
and energy for long term integration with respect to time 
of tha governing time- dependent equations . 

2. 6. 4. 2 Finite Difference Equations 

The rar-field and near-field versions of the rigid-lid 
model use the same set of finite difference equations; although 
the initial conditions and boundary conditions are cuite 
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different. Equations (2. 9) -(2 13) for continuity, hori~ 
zontal momentum conservation, and conservation of energy are 
approximated by finite difference equations by using a for- 
ward time, central space scheme (so-called FTCS) , with DuFort- 
Frankel (1953) differencing performed on the diffusion terms. 
The finite difference equation for u-momentum may be written 
as for example 

— 2 ;— (Inertia)*^ + (Coriolis)’^ + (Pressure)*' (2.30) 

+(Viscous)’''’'‘^^'***'^ 

wh'"*e the DuFort-Frankel differencing is expressed as 

_3fu . u?I+I,J.K)-l- U°(I-X,J,K) - u"~^^(I,K,K)-u"~^g,J,K) (2.31) 

(.V.)’ 

The effect of modified DuFort-Frankel differencing is to relax 
somewhat the diffusive stability criterion, Sengupta (1974). 

The surface pressure equation (2.17) may be solved by 
iteration, at each time step, by using successive over- 
relaxation or by the modified SOR technique (Liebmann Method) . 

Fig. 2.4 shows the flow chart for the steps involved 
in propagating the numerical solution of the system of govern- 
ing equations for the rigid-lid model. These steps may be 
elaborated as follows : 

1. The problem is set up as an initial value problem. 

The values of u,v,P,p,il and T are specified initially for 
time level n. 

2. Using the known values of the system variables at 
time level n the forcing function, Tr(a,3) is evaluated at 
all half-grid nodes. 

3. Surface pressure, P , at the "rigid lid" is eval- 
uated by iteration (Liebmann method) at all half-grid nodes 
using the Poisson equation. 
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4. Preisure gradients computed from the Poisson 
equation are used to compute u and v by the horizontal mom- 
entum equations for the current time level n + 1. The 
hydrostatic equation (2.12) is used in the process to 
obtain the three-dimensional pressure field. 

5. The continuity equation is then used to compute 
equivalent vertical velocity, , at time level n + 1 from 
the known values of u and v at time level n + 1. The values 
of are obtained by integrating the continuity equation 
from Y-0 to Y*l. 


6. The actual vertical velocity, w, is then computed 
at time level n + 1 by using equation (2.16). 

7. The energy equation is then used to compute the 
values of T at time level n + 1 from the known values of 
u,v and at time level n + 1 and T at time level n. 

8. The solution domain is then checked for static 
stability. If there is cooler water on top of the lower 
warmer water, i.e., if 

%4l| J+^ for unstable (2.32) 

conditions infinite mixing is invoked. 

9. The density, P , is then computed from the quation 
of state knowing T at time level n + 1. 

These nine steps are then repeated to propa gate the numer- 
ical solution to time levels n + 2, n + 3, etc. 

2.6.5 Stability Criteria 

Since it is not possible to make a strict stability 


analysis for the system of governing equations under con- 
sideration, the one-dimensional iur^;erpequation is used for 

stability analysis. This is relevant since the Burgers 
equation contains an unsteady term, a convective term and a 
diffusion term. The stability criteria for'iuraers equation 


+ r 

3t ^ ^x 


3u 

3x 




(2.33) 


as discussed by Roache (1972) are as follows 
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CONVECTIVE : 


(2.34) 


C 


X 


.\x 


I 


DIFFUSIVE: 0 i 

^ (-Tx)^ 


(2.35) 


Th« Burgers equaclon represents the one-dimensional form of 
incompressible fluid motion. The convective stability 
criterion may be interpreted to require that no fluid 
particle can move the distance of one spatial grid space 
in one time step. Likewise, the diffusive stability 
criterion may be interpreted to require that momentxim 
cannot diffuse to half the distance of one spatial grid 
space in one time step in forward and backward directions . 
Thus, for the numerical solution to be stable, the time 
step must be small enough to give sufficient time for the 
physical processes to develop at each grid node. 

The stability criteria may be extended to the three- 
' 'mensional equations as follows: 

CONVECTIVE: ^ + C^^At*^ • 1 (2.3b) 


DIFFUSIVE: ut 0 + K ■) + D„ .'t • 'n (2.37) 

^ (3x)" y CJ)" “(!?)“ 

For the application of these criteria to the present problem, 

C.^» C,T> C.T ^7 interpreted as the maximum values 
X y z 

of u, V, and w in the domain; and D^^, D^, and D^ may be 
interpreted as the kinematic eddy diffusivities in the 
X, y, and z directions, respectively. 

2.6.6 Higher Order Terms 

Higher order terms resulting from the transformation 
of the horizontal diffusion terms (i.e. second order 
derivatives in‘^ and 3 ) , from the Cartesian coordinate 
system (x, y, z) to the vertically stretched coordinate 


58 


8yst«m (ot ii Y ) have been neglected. This has been done since 
the magnitude of the vertical diffusion terms are several 
orders of magnitude larger than the horizontal diffusion 
terms, Sengupta and Lick (1974). Appendix A of this volume 
presents the details of this transformation. 

2.7 Sample Results 

In this section sample results using the rigid-lid 
model for near-field and far-fleld applications will be 
presented. 

2.7.1 Near-Field, Cutler Ridge Site (open basin) 

The region of Influence of the Cutler Ridge plume has 
been approximated to a rectangular domain as shown in 
Fig. 2.4 which is open on three lateral boundaries. The 
grid system for the rigid-lid near-field model of the 
Cutler Ridge site is shown in Fig. 2.5. The rectangular 
domain of solution extends 425 meters in lateral extent and 
525 meters in longitudinal extent. The discharge is taken 
as 25 meters wide. The nume-**ical grid system has 18 and 
22 nodes across and along Che axis of Che jet, respectively. 
There are 5 nodes in the vertical direction. 

In order to understand the physical processes involved 
and to investigate Che numerical behavior of the model 
several simplified cases were executed before the final 
calibiation and verification run for the data bases obtained 
for April 15, 1975. Volume II presents the details of this 
computer run. Table 2-1 shows different cases studied 
together with important features for each case. Lee and 
Sengupta (1977) . 

The data base for the final calibration and verif- 
ication run is obtained from the field experiments conducted 
at the Cutler Ridge site on April 15, 1975. The initial 
temperature conditions are taken from the morning IR data 
(0911-0912 EST) shown in Fig. 2.6 and ground truth data on 
April 15, 1975. The computations were continued for 2^ 
hours. Figs. 2.7 to 2.10 show velocity distribution at Che 
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surface and at different vertical layers. The interaction 
between the wind driven current and the plume due to current 
coming from the south can be seen in these figures . The 
fluid is exiting through the east and western boundaries 
only. The reduction of the effects of wind ar>d current in 
the lower vertical layers can also be seen due to bottom 
shear. Fig. 2.11 shows isotherms predicted by the model 
along with IR data (11:45 - 11:55 EST) on April 15, 1975. 

As can be seen, there is reasonably good agreement between 
IR data and model predicted results. Fig. 2.12 shows 
temperature decay along J at I«ll (i.e. close tc the center- 
line) predicted by the model along with IR data. There is 
good agreement between IR data and model results. The iso- 
therms along I and J sections are shown in Figs. 2.13 and 

Stratification near the discharge with 
isothermal conditions away from the plume can be seen in 
Fig. 2.13. 

2.7.2 Far-Field, Biscayne Bay (closed basin with ocean 
efflux) 

Ignoring the Cutler Ridge site thermal discharge, the 
rigid- lid far-field model was applied to Biscayne Bay to 
investigate the naturally occuring circulation and far-field 
temperature distribution. 

Applying the rigid-lid far-field model to Biscayne Bay, 
solutions have been obtained in three stages. First, a 
closed basin approximation was made and wind driven cir- 
culation was obtained. In the second stage circulation in 
the Bay was obtained with an ocean efflux specified; the 
simultaneous effects of wind and ocean efflux were also 
investigated. The third and final stage obtained the temp- 
erature field for various ambient condition*’ and parameters . 
The results were compared with airborne thermal scanner 
IR data t;> calibrate the model. This procedure was followed 
in order to check the performance of the model for as wide 
a range of environmental situations as possible. The hor- 
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izontal grid system for Biscayne Bay is shown in Fig. 2. 15 
(I-ll, J»34) . Five vertical layers were used in the model. 

Lee and Veziroglu (1975) present results for this 
application for the three stages noted. Only results for 
the third and final stage will be presented and discussed now. 

A 10cm/ sec velocity for both the incoming and oxitgoing 
tidal phases was assumed. The program is executed with an 
incoming tide, and then the ocean-bay velocity pattern is 
gradually reversed to obtain a lOcm/sec outgoing velocity. 
Although the details of phase relations, level changes and 
time dependence cannot be precisely modeled by the rigid- 
lid model (c.f. section 2.1 of Volume I), the results should 
give a meaningful equalitative picture of the circulation. 

Fig. 2, 16 .hows the surface velocities with an in- 
coming tide. The major portion of the tidal mass influx 
travels into the South Bay while the flow towards the 
closed northern region is minimal. Fig. 2. 17 shows the 
velocities at a depth of 1 meter with an incoming tide. 

The flow is unidirectional at most points . The incoming 
flow was gradually reversed. Fig. 2. 18 shows the surface 
velocities at an intermediate stage. The currents have re- 
versed in some places but not in others. In Fig. 2. 19, the 
currents at a depth of two meters are almost completely 
reversed. 

Fig. 2.20 shows that the bulk of the outflow comes 
from south bay. At Cutler Ridge the current is now from 
west to east, therefore, with outgoing tide the plimae is 
expected to turn towards the east. 

The effects of wind and tide on the general circulation 
are shown in Fig. 2.21. The effect of the southeast wind is 
to turn the current vectors toward the northern part of 
the bay. The velocities in the south bay are still essentially 
southward but are of a smaller magnitude due to the south- 
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east wind. 

The temperature solutions for various combinations of 
parameters have been obtained for comparison to the April 
15, 1975 IR data base. Table 2-2 shows a list of some of the 
cases for which solutions were obtained. Fig. 2.22 shows 
good agreement between model rerv.lts and the IR data base 
for April 15, 1975. Fig. 2.23 shows vertical section J»7 ill- 
ustrating the isotherms in a transect. As can be seen 
vertical temperature variation is relatively small owing 
to the shallowness of the bay and the turbulent mixing 
processes . 

2.7.3 Near-Field and Far-Field, Lake Belews 

For ease of mathematical modelling, the total path 
of water circulation is divided into two regions , namely 
the mixing pond (near-field) and the main lake (far-field) . 
These two regions are treated as disconnected regions. 

Fig. 2,24. shows the grid system for the mixing pond (1*29, 

J»13) . Fig. 2.25 shows the grid system for the main lake 
(I»29, J“13) . Lee and Sengupta (1977). 

The mixing pond receives hot water from the power 
plant, mixes it with cooler water and then discharges it 
into the main lake through the connecting canal. The main 
lake receives hot discharges fj/om the connecting canal, 
cools it and from there it goes into the power plant condenseis. 

The primary difference in the parameters in the mixing 
pond and the main lake is the fact that the mixing pond is 
well mixed while the main lake shows thermal stratification. 
This in turn means that the vertical eddy diffusivity 
in the mixing pond is constant over the entire depth whereas 
in the main lake vertical diffusivity decreases with depth. 

Table 2-3 gives th-i list of computer runs for Lake 
Belews site. Fig. 2.26 shows the surface velocity pattern 
in the mixing pond for the August 23, 1974 data base; and 
Fig. 2.27 shows the velocity distrubtion at four meters 
depth in the mixing pond. It can be seen that at the 
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surface the current vectors are in the direction of the wind 
to right of the connecting canal. Whereas at four meters 
depth the circulation to the right of the canal is reversed 
from that on the surface. Figs. 2,28 and 2,29 show surface 
isotherms and isotherms at four meters depth in the mixing 
pond for the August 23, B74 data base. Fig. 2.30 shows 
a comparison of model predicted surface isotherms and IR 
data base for May 19, 1976. The general agreement is 
reasonably good. 
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Fig. 2 - 10 ' Velocity Distribution for April 15, 1975 
at K=4 for Cutler-Ridge Site (Case NO, 8) 
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Fig. 2-15 Horizontal Grid System for Biscayne Bay 
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III. 


The Free-Surface Model 

3 . 1 Brief Description of Past Experience with Free- 
Surface Formulations , Relevant Advantages and 
Disadvantages 

One of the first three-dimensional models was by Freeman 
et al (1972) . This model was essentially a free surface 
formulation. Haq and Lick (1974) used a free surface 
model to study the time-dependent flow in large lakes with 
application to Lake Erie. They used a vertical stretching 
to convert a variable depth basin to constant depth, thereby 
permitting a constant vertical grid size to be applied 
everywhere in the domain. Irregular shorelines could be 
easily included without modification of the computer 
program. Modified versions of the free-surface model have 
been used, by the Thermal Pollution research team at the 
University of Miami, at a number of sites with satisfactory 
results. Lee and Sengupta (1976). 

The major advantage associated with the free-surface 
model is its ability to predict surface heights everywhere 
in the domain. Thus, for example, real tidal conditions 
can be accounted for by this model, and, hence, verification 
of the model can be made with regard to comparison with 
existing tide data bases. 

The major disadvantage associated with the free-sur- 
face model is its inherently small time step as determined 
by the Courant-Freidrichs-Lewy Condition^ Roache (1972), 
Richtmyer and Morton (1967) , which is based upon external 
gravity waves ( or so - called surface gravity waves) . 
However, for water bodies for which vertical diffusion 
determines the maximum allowable time step,, as. would be 
the case for shallow water basins , there is no real time 
step disadvantage in using the free-surface model. 

3 . 2 Assumptions and Approximations 

The system of governing equations (see next section 3.3) 
for the fluid flow invoke several simplifying assumptions and 
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approximations in the interest of saving computational time 
without losing significant accuracy. The following assumptions 
and approximations have been employed. 

3.2.1 The Boussinesq Approximation 

The effect of density variations on the Inertial and 
diffusion terms in the governing conservation equations is 
neglected. Density variation is retained in the bouyancy 
terms in the equations of motion. The effect of bouyancy 
is thereby accounted for by allowing density variations 
which produce horizontal pressure gradients which influence 
the fluid motion through the horizontal momentum equations . 

3.2.2 The Hydrostatic Approximation 

The hydrostatic approximation involves neglecting the 
vertical convection and diffusion terms in the vertical 
momentum equation. This approximation implies that the 
vertical fluid acceleration,^ , is negligible. 

3.2.3 Constant Eddy Transport Coefficient s 

Turbulence modeling is very complex and has an ex- 
tensive body of literature of its own. Turbulent closure 
has been obtained in this model by using constant eddy 
transport coefficients, although the horizontal eddy trans- 
port coefficient is orders of magnitude larger than the 
vertical eddy transport coefficient, being due to' the much 
larger horizontal scale length, L, in c onparison with the 
vertical scale length, H. 

3.2.4 Variation of Surface Wind Stresses 

The variation of the wind produced surface shear 
stresses with respect to x and y, |J^^ and |- T7? , are con- 
sidered negligible for the horizontal length scales of the 
water bodies studied. However, if the physical dimensions 
of the water body are so large as to require including 
variation of the wind stresses with respect to x and y, 
then the computer programs can be quite easily modified by 
replacing t andx with matrices t (I,J) and t (I,J) 

Where the indices I and J refer to the location of a grid 
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point with respect to the x, y plane. 

3.2.5 Velocity S l ip Conditions 

The free-surfac;,, far-rield model uses velocity slip 
conditons at the lateral solid boundaries , although no-slip 
is used at the bottom boundary. The assumption of slip 
conditions is necessary for the free-sxirface model to allow 
for surface height variations at the solid boundaries, 

Freeman et al (1972) , Lae at al Cl 97 6) Numerically it has 

been seen that lateral boundary layers are smaller than the 
relatively large grid spacing used, specifically for the 
Biscayn- 3ay site. Estimates by Sengupta and Lick (1974) have 
indicated that the sidewall boundary layers are thin for 
similar situations, and do not extend as far as the nearest 
interior node. Note, that the free-surface, far-field 
model uses the same velocity slip conditions used by Freeman 
et al (1972) . 

C*) Specifically , at x-bomdaries «o , and at y-boundaries 

3 (Hu) - 0 . “ 

ae 

3.3 Governing Equations 

The set of equations governing the behavior of the 
fluid flow are those expressing the conservation of mass, 
moment\am, and energy in turbulent flow, and an equation of 
state. 

3.3.1 Cartesian Coordinate Representation (x.y,z) 

The Cartesian coordinate system is used with the 
z-coordinate in the downward vertical direction as shown 
in Fig. 3.1, i.e. a so-called "left-handed" coordinate system. 
In order to keep the generalized nature of the model, all 
the significant terms in the respective conservation equations 
are retained. Included are the effects of bouyancy, inertia, 
Coriolis , density and turbulent mixing . Wind shear and heat 
flux at the surface are also considered. 

(0 Note : ~ ’ »o and ^ =*o yields "verv~close results. 

3a. .33 . ..r: — - 
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The following system of non-linear partial differential 
equations, written in Cartesian coordinates, describes the 
three-dimensional, unsteady fluid flow where the variables 
are in dimensional form, 

Continuity Equation 

3u + ^ + ^ ■ 0 (3.1) 

^ 3z 3z 


Momentxjm Equation 


3u + . 3u + 


It 


u- 


3x 
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3u 
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if 
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3y 
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+ Kjj l!v + i!v 

3y 


3y 


p 3z ^ 


3x 


3y 


(3.3) 


+ K„ 3^w 

V y 

3z 


(3.4) 


Energy Equation 

3T + 3T + 9T + 3T = 


3t 


u 


3x 


V 


3y 


w 


9z 


Bg ifr + Eg sfr + ifi 

9x 9y 9z 


(3.5) 


Equation of State 

p = p(T) (3.6) 
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Continuity Equation 


3H + 3iHu) + 3 (Hv) + H 3S2 - .0 
3t 3a SS 3a 


(3.9) 


Horizontal Momentum Equations 

3 (Hu) + 3 (Huu) + 3_(Hw) + H 3 (ufl) 

3t 3a '3^ 30 

■ « C- 5 + 8 (“I? - 15 ) ^ fv] 

It ^ ■*■ p ^ If ^ 

(3.10) 


g_(Hy)- + 3 (Huv) 3(Hw) . „ 3 (vO) 

at 3a 33 30 

- « C-i (||) + 8 (oH - If) 


+ Kn[b («lf) +% 


^ /Tj a V \ 1 1 a 

a3 3 3'' p HTo" 


PK 


V 3 
(3.11) 


3 0''J 


Energy Equation 


a (HT) . 
3t 


a (HuT) , 3 (HvT) 
3a 3 3 


+ H 


3 (ni; 


30 


= B. 


H 


{h li)] 


®H Cl B 




H3 0 3 

(3.12) 


i>] 


Equation (3.4) for conservation of w-component of tnomentum is 
replaced in the free-surface model formulation by applying 
the hydrostatic approximation (see section 3.3.2) as follows: 
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(3.13) 


Hydrostatic Equation 

a 

P(a) - P (0 -0) +gH / (a) da 

0-0 

Equation of State, p*p(t) is given for fresh and salt water 
as follows : 

Saltwater; p(T) -.1.029431 - .000020T - .0000048T^: .. 

Cfor a salinity of 38 parts per thousand) (3.14) 

Fresh Water :p(T) - 1.000428 - .000019T - .0000046X^^3 


Instead of using equation (3 9), following the work by 
Freeman et al (1972) , two integrated forms of the continuity 
equation are used as follows : 

Surface Height Equation 


3_H 
3 t. 


- - + S.(.Hv) ^ 

^ Caa ^38 J 


(3.16) 


0»O 


Equivalent Vertical Velocity ('ina,S,0 system^ 

_ _1 r 3 (Hu) . 3 (Hv) . 

“ I '■ ^ 3l ^ 

0-0 

1 


+ 0 ; ^ 3iHuI do 


H 0=a 


36 


Actual Vertical Velocity (in x,y,z system) 


w » HJ2 +0 ^ + ( 0 - 1 ) ^ 


(3.17) 


(3.18) 


where, ^ 


3h + . 3h + 3h 


3 t 


‘3 a 


Vi 


36 


dn 

It 


3n + u3n +-3n 
3 1 3a ^38 


and, 


w 


dz 

It 




d0 

Tt 


(1) Note: These two integratt'd forms of equation (3.9) 

consider ^-0 at '7=0, and f^-0 at 7=1. However, a special 
case of mass influx at the bottom boundary, where at 
0=1 is given in Volume III for the sample problem for 
Hutchinson Island site. 


The symbols in equations (3.1)- (3.6) are defined in the 
list of symbols for the free-surface model. 

3.3.2 Vertical Stretched Coordinate Representation (a,3,c) 

One major difficulty in the treatment of the free-sur- 
face model is at the free surface boundary. The boundary con- 
ditions can be specified, but the position of the free sur- 
face is irregular and time-dependent making it very diff- 
icult to apply any grid system at this boundary for numerical 
solution. The approach used in the model formulation is to 
follow a vertical stretching transformation suggested by 
Phillips (1957) and used successfully in studies by Freeman et 
al (1972) . Using this transformation, the free surface be- 
comes a fixed flat surface and the variable depth u ^,om be- 
comes a flat bottom boundary. This method allows easy adapt- 
ation to various bottom topographies, an important require- 
ment for any general model. In addition, constant vertical 
grid size can be used throughout the domain. 

The transformation of the vertical coordinate for the 
free-surface model is obtained by letting 

a = X 

6 = y (3,7) 

and a = Z(x,y, z, t) = z+ n(x,y, t) 

H(x,y, t) H(x,y, t) (3.8) 

'where the symbols are given in the list of symbols for the 
free-surface model. Fig. 3.2 shows the (a, 8, a) coordinate 
system. Note, that the value of a ranges monotonically from 
zero at the free surface to unity at the bottom boundary. 

By substituting transformations (3.7) and (3.8) into 
equations (3.1) - (3.6) the free-surface model governing 
equations (in dimensional form) in the (a,3,cf) coordinate system 
are expressed in what follows. 
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3.4 Horizontally Stretched Equations 


It is desirable to obtain a more detailed description of the 
flow near the discharge point while larger grid size may be used 
in the further points to save computation time. A horizontal 
stretching of the coordinate system (a, 3, a) is used here to 
create a more efficient use of the grid points by letting 

a - a + Sinh ^ C 2 (X-d) ] (3.20) 

3 » b + C 2 Sinh (Y-e)^ 


where the various symbols are defined in the list of symbols 
for the free-surface model. Fig. 3.3 shows the X, Y,<T coordinate 
lystem, and Fig. 3.4 shows the resulting u, 3 ,0 , coordinate oystem. 
^Ippendix B presents the details of this coordinate transformation 
and the resulting equations in the coordinate system. Also 
comparison of sinh stretching with tanget stretching, used by 
Waldrop and Farmer (1973) is given in Appendix B. 


Now, after defining the following derivatives necessary for 
making this transformation, the horizontally stretched, free- 
surface model equations will be presented. 


dX 

3a 


Y’ 


d3 ' 


d^X 


t 



The transformed free-surface governing equations in the 
(X, Y> a) coordinate system are as follows: 


Continuity 

3H + 3 (Hu) + 3 (Hv) + „ 3 

3t ^ 3X ^ 3Y ^ fa 


(3.21) 


Horizontal Momentum Equations 


3 (Hu) _L v' 3 (Huu) 3 (Huv) , „ 3 (ufl) 

“TT” ^ ^ ^ ~ 1 ~ ^ TY" ^ ^ “33~ 


3u 

3X 


, Ir 1 9 

+ D '••H 5^ 


9H 

9u 

9: 

■T. : 
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. a 

9Y 

9Y 

( 

9u 

9a 



9 







The hydrostatic equation and the two equations of state (for salt 
water and fresh water, respectively) are still given by equations 
(3.13), (3.14) and (3.15). 

Surface Height Equation 


3 t 


1 

, t ligui ^ y’ lOW 3 

a-0 9Y 


C3,25) 


For, the free-surface model application to the sites where the 
velocity on the bottom of the basin is not zero, such as a 
submerged discharge the vertical velocity on the bottom 

no longer is zero, and then the integration of the continuity 
equation to get the equation for surface height can not eliminate 
the vertical convection term. An integration constant re- 
presenting vertical velocity at the bottom, is added to the 
equation. Thus, from the relationship of w and we get: 

n, = 1 (W, - u.X 3h - V Y-* 3h ) (3.26) 

^ H ^ ^ ^ 3Y 


where subscript b denotes the fluid property at the bottom of 
the basin. By integrating the continuity equation (3.21) with 
respect to cr from the free surface (c-o) lo the bottom (a=l) 
we get: 

3H _ r r X 3 (Hu) +Y 3 (Hv) i da 
pt “ “ 3X ~TT~-^ 

a =0 

- <«b - “b H -''b Iy > 

^^^NOTE: This is for the case =0 at a=l. 
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Equivalent Vertical Velocity 


(3) 


SI 


..1 do 

a=* 0 

+ 0 1 
H ^ 


C x' i%l + 1^1 da 

q = d 


(3.28) 


For a submerged discharge at the bottom, substituting equation 
(3.27) into the continuity equation (3.21), and then integrating 
from the free- surface to we get for the equivalent vertical 

velocity : a 


SI 


a*0 


H /. 


3 (Hu) 
ax ‘ 

3 (Hv) . 
aY ] 

da 

3 (Hu) . „• 

ax ^ 

3 (Hv) , 
aY 

da 


3h 


+ H »b - “b O 


V v' 3 h \ 
b * §T-' 


(3.29) 


3.5 Boundary Conditions (closed basins and open basins) 

The nature of the system of governing equations requires 
initial and boundary conditions to be specified. The boundary 
conditions for both near-field (open basin) and far-field (closed 
basin) versions of the free-surface model are presented in this 
section. The initial conditions will be presented in the next 
section. 

3.5.1 Combined Near-Field and Far-Field (open basin) 

The. free-surface model with horizontal stretching has been 
applied to a submerged thermal discharge into the ocean at a 
coastal site (open basin) . Tsai (1977) . The boundary conditions 
are in summary : 

*'^^NOTE: Again this is for the case S>= 0 at 1. 
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I 


At the 

Surface, a » 0 

At lateral solid boundaries 



on x-boundaries 



^ » 0 


n 0 



3U 
3 a 


u » 0 

V » 0 



3v « - 
3a 

^zy 

3T , 3T , X-' 

^ 3X ^ H 

3n 

3X 

3T 

3a 

3T - - 
3a 


y' a 3 H 3 T 
H 3X 3a 

a 0 


At the 

bottom, a *1 

on y-boundaries 



* 0 

(except at sx±merged discharge) 

Q. f 0 


( 

Ij 

u - 0 

(except at submerged discharge) 

u* 0 


i 

1 

V a 0 

(except at submerged discharge) 

V - 0 


i 

0 

3a 

(except at submerged discharge) 

aT»y’ 3T +t’ 3n 
Iy ^ 3Y H 3Y 

3T 

3a 

t 


At lateral open boundaries 
on x-boundaries 


Q 


3u 

3 X~ 3X"^ a ax 3 q 

f' 

- X £ 3H la =0 
H 3X 3a 


£ ^ 3T =» 0 
H 3Y 3a 


on y-boundaries 
0 


^ §L3i 

3Y ^ 3Y H 3Y 3a 

- Y 'a d_E 3u = 0 
iT 3Y 3a 


(4) HOTE: H is depth contour for free-surface model. 
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4 


3v 

3X 

Y ' 9v 
^ 3l^ 

+ 

X- 3n 
H 3X 

3a 

3v 

r? 


3v , 
SY'^ 

Y 

5 

9ri 

3Y 

3a 


Y ■ '0 

3H 

3v ■ 

0 


1 

a 3H 

3v 

as 

0 



3a 





3a 




0(')2|S II + h(x')2»^ 




wz 


+ HX'- |T;|.o 




+ Hy'||].o 


at Discharge 
Velocity ^ 
Temperature ^ 
Density I 


Specified 


3.5.2 Far-Field (closed basin) 

The free-surface , far-field has been applied to the Biscayne 
Bay, Lee et al (1976). The bo'undary conditions are in sxammary: 

At the surface, a= 0 

= 0 


(5) At lateral solid bou ndaries 
on x-boundaries 




H 




LI 

3a 




n i>0 
u =0 


3CHv) ^ « 
3a ^ 


= 0 
3 a 


P V 

At the bottom, a = 1 
= 0 
u = 0 

V = 0 


o n y-boundaries 


7^ 0 

9 (Hu) ^ . 
3 B 

V = 0 


3T _ 


3 R 


= 0 


0 
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At lateral-open boundaries (current velocity specified at 
inlet) 

OK x-boundaries on v-boundaries 

ft 0 ft f 0 


( 6 ) 


u ■ u (t) or 


3u 

da 


V “ 0 



u » 0 

V ■ V (t) or - 0 


At lateral open botindaries (tide height specified at inlet) 
on x-boundaries. on y-boundaries 


ft # 

0 or^^i f 0 

ft 0 or^^^ 

ft 0 

(U) 

n ■ nCt) 

u - 0 

u - 0 

3u 

3a 

0 1^-0 
3a 

-n(t) 


V » 

0 V » 0 



3 a 

0 H-0 


5-1 . 0 

9S 

3.6 

Initial Conditions 




The initial conditions are specified by using the corrected 
doming IR data base for temperature; as will be illustrated in 
the sample problems in Volume III and zero velocity everywhere in 
the domain (u * v * ft » 0) , since it is quite difficult to 
obtain ground truth current measurements for the entire domain 
for the kind of grid size resolution that would be required. 
Although an initial free-surface n*nCx,y, t»o) ,can be specified 


NOTE: ft-sois used in program to save computational timej whicn 

is a good approximation, since it has been learned that 

NOTE: At an inlet uCt). or v(t) must be specified, and at an 
outlet ^ or ^ ^ay be specified. 

3a 3 3 


1G9 


from existing tide data bases, compatibility between the sur- 
face heights and the velocities requires starting the computations 
with a flat surface, n(x,y , t-0)-0, initially. 


^^^NOTE: n*<»is used in program to save computational time, 

which is a good approximation, since it has been learned thatJ?*o- 

3 3v 

'' ''NOTE: At an inlet n-n(t) and either are specified; 

and at an outlet 3 u or- a v specified. 

STcT 515’ 
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3 . 7 Method of Solution 


3.7.1 General 


It is obvious that closed-form analytical solution of the 
system of governing equations (3.9) - (3.18) for the free-surface 
model in the(a ^3, cr) coordinate system, and for equations (3.21)- 
(3.29) in the (X,Y,a; coordinate system, is impossible to get. 

The set of equations consists of coupled, unsteady, three-di- 
mensional, nonlinear partial differential equations. Therefore, 
the finite difference method is used to ob tain numerical solutions . 

A three-dimensional grid system is established with respect 
to the^o»,^3; cr) coordinate system for the vertically stretched 
free-fei^f^ane model equations, and a three-dimensional grid sys- 
tem is, e‘st‘ablished with respect to the(X,Y,a) coordinate system 
for the vertically and horizontally stretched free-surface model 
equations. The governing equations are then solved over finite 
time steps which are carefully selected to obey numerical stab- 
ility criteria. This will be discussed in detail in a following 
subsection on stability criteria (3.7.5). 

In general, several methods are available for integrating 
over time the governing equations for incompressible fluid flow 
as discussed by Roache (1972) . The two most common techniques 
for integrating time-dependent partial differential equations 
are the implicit and explicit finite difference methods . The 
implicit method involves the solution of a set of simultaneous 
equations which are obtained by writing the spatial derivatives 
in terms of the respective unknown quantities at the current 
tim-B level n + 1, knowing the values of the remaining quantities 
of the set (u, v, , h, P, p, T) at the two previous time levels 
n and n-1. For the one-dimensional case, implicit methods are 
convenient because the set of simultaneous equations is tri- 
diagonal, Richtmyer. and Morton (1967), and, hence, a direct 
matrix inversion method of solution is used. However, in the 
case of a three-dimensional model, the implicit method becomes 
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too time constimiiig , since the simultaneous equations must be 
solved at each time step by an iterative technique. Thus, 
although the advantage of implicit methods is that they allow 
larger time steps, for the three-dimensional case the iteration •’ 
time for each time step more than offsets the inherently larger 
time step. Furthermore, alternating direction-implicit (ADI) 
methods may be used to obtain tri-diagonal matrices even for 
multidimensional equations, however, for irregular boundaries 
the ADI methods are impractical. 

Therefore, the explicit finite difference method is used 
for numerical solution of the free-surface model. The solution 
to a particular partial differential equation is propagated from 
point to point on the numerical grid system. The current time 
level value., n + 1, of a particular system variable (u,v^H, P,p, T) 

is computed in general from known values of the corresponding 
system variables at the two previous time levels n and n-1. 

Thus, this is an explicit scheme. 

The mathematical model is an initial-value, boundary- 
value problem and, hence, requires specification of both initial 
conditions and boundary conditions (see sections 3.5 and 3.6). 

3.7.2 Computational Grid System 

The free-surface model does not use the staggered grid 
(or mesh) system as used in the rigid- lid model (section 2.6.2). 
Instead, the full grid system is used for defining the system 
variables u,v,n,H,P,p,T at the integral nodes (I, J. K) . The 
rigid-lid model uses the half-grid system for better meshing of 
the solution of the Poisson equation for surface (or lid) 
pressure, P , with the horizontal velocity components u and v. 

This is not considered necessary for the free-surface model, since 
the pressure field P (I, J, K) is computed from H (I, J) which 
is computed at integral nodes from u (I, J, K) and v (I, J, K) . 
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3.7.3 MAR Ntimbarlng Systain 


Since the free- surface model does not use the half-grid 
system, only MAR (I, J) is used for distinguishing between 
spatial finite differences in the interior: on the boundary, and 
outside the domain of solution. The MAR numbering system for 
both the far-field (closed basin) and for the combined far-field 
and near-field (horizontal stretching applied to an open basin) 
versions of the free-surface model will be clearly specified in 
the sample problems in Volume HI. Note, that for the far field 
model MAR (I, J) "6 and MAR (I, J) » 8 boundary corners are 
treated as interior points (MAR(I, J)-ll) 

3.7.4 Finite Difference Schemes 

3. 7.4.1 Approximation of Spatial and Temporal Partial Deriv- 
atives ; Conservative Form 

The spatial derivatives are central differenced in the in- 
terior, for example: 

a_u - ud-M.J.K) -u(I-l.J . Kl 

9 a i4a " (3.30) 

and, 

i u(I+l,J,K) +uq-l.J.K) -2u(I.J.K) 

3a^ r4a)z (3.31) 

At the boundaries , three-point single sided schemes are used 
by fitting a parabola through three points (the boundary point 
and the next two coincident interior points). Thus, for example, 
at the left a -boundary; 

Lu = 4u(I+l,J.K) -3u(I.J.K) -u(I+2.J,K) 

9a 2 a a (3.32) 


113 


u(I.J.K) +u(I+2.J.K^ -2u(I+L,J.K) 


and, at the right a -boundary 

iii - 3u(I.J.K) + UCI-2.J.K) -4u(I-l.J.K 


and, u(I.J.K) +u(I-2.J,K) -2ug-l.J.K) 


Kota, that the spatial finite difference approximations (3.30) 
(3.35) are on order of accuracy of ct )?Crandell (1965). 


The temporal derivatives can be expressed in two forms, 


first 


3u « uj 
3 1 


- uu.. 


(3.36) 


for forx^ard differencing in time, which is on the order of 
accuracy of At and 

n+1 n- 1 

3u » u(I.J>K) - u(I,J,K) 


for central differencing in time, which is on the order of 

2 

accuracy (At) . 







» 

1 

i 

The finite differences in both space and time in the model 
are expressed in the fi\ll conservation forms following Arakawa 
(1966), for example: 


This is done to avoid possible "leaking" of mass, momentvim, and 
energy for long term integration with respect to time of the 
governing time-dependent equations , 

3. 7.2.4 Finite Difference Equations 

The full set of finite difference equations for the far- 
field and combined far-field and near-field (horizontal stret- 
ching) versions of the free-surface model will now be discussed. . 

3. 7. 4. 2.1 Far-Field 

The two integrated forms of the continuity equation for the 
surface height H, (3.16) and for the equivalent vertical velo- 
city J?. (3.17), are integrated over depth by applying Simpson's 
rule. The time derivative in the surface height equation is 
initially replaced by a forward difference in time, and there- 
after a central difference in time is used. 

The numerical method used for solving the horizontal mo- 
mentvrai equations for u and v is an explicit finite difference 
scheme for which a forward difference in time and central 
differencing in space (so-called FTCS) is used. The horizontal 
diffusion terms are differenced at n-1, i.e., two time steps 
back from the currently computed time level, n+1. The vertical 
diffusion terms are differenced using the DuFort-Frankel scheme, 
Roache (1972) . The V!tCS method is used througout for solving 
the energy equation without Du-Fort-Frankel differencing of 


3 (Hu ) ^ 
3a ” 


(Hu) 


I-l, J,K 


(3.38) 
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the vertical diffusion term., and has produced no ntjmerical . . 
instability problems. 

3. 7. 4. 2. 2. Combined Far-Field and Near-Field (Horizontal 
Stretching) 

Spatial integration of equations (3.25) and (3.28), or 
of (3.27) and (3.29), have been performed by applying the Trap- 
ezoidal rule. Again, FTCS is used for solving the u and v 
momentimi equations with the horizontal diffusion terms evalu- 
ated at n-1 and the vertical diffusion terms DuFort-Frankel 
differenced. The energy equation, also uses FTCS with the 
same differencing of the respective diffusion terms as done 
for the horizontal momentum equations. 

3. 7. 4. 2. 3 Flow Chart 

Now we see the flow chart for the steps involved in pro- 
pagating the ntmierical solution of the system of governing 
equations for the free- surface model. These steps may be elab- 
orated as follows: 

1. The problem is set up as an initial-value problem. 

The values of u,v,n,H,P,p, and T are specified initially for 
time level n. 

2. The surface height equation is then used to compute H 
at time level n+1 from the known values of u and v at time level 
n. P is computed at time level n+1. 

3. The horizontal momentum equations are used to compute 
u and V at time level n+1 from the known values of u,v,^2,p,p 
at time level n. 

4. The equivalent vertical velocity equation is used to 
compute at time level n+1 from the known values of u,v 

and H at time level n+1. 
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5, The energy equation is used to compute T at time level n+1 
from the known values of u,v, and H at time level n+1. 

'6. The actual vertical velocity, w, is then computed at time level 
n+1 from u, v and H at time level n+1, 

7. The density, p , is computed from the equation of state know- 
ing T at time level n+1. 

These seven steps are then repeated to propagate the numerical 
solution to time levels n+2, n+3, etc, 

3.7.5 Stability Criteria 

The one -dimensional Burgers equation is used in the stab- 
ility analysis by a heuristic extension into three-dimens ion s- 
This method follows Roache C1972) , A strict stability analysis 
for the system of governing equations under consideration is not 
possible, The stability criteria may be extended to the three- 
dimensional equations as follows; 

CONVECTIVE: (|^) + C^(^) + < 1 (3.39) 

DIFFUSIVE : (Ax)^ + D„(Ay) + D^(Az)^< % (3.40) 

X y z 


For the application of these criteria to the present problem, C , 

C , and C may be interpreted as the maximum values of u,v, and w 
y z ^ 

in the domain and D , D , and D may be interpreted as the kine- 

X y z 

matic eddy diffusivities in the x,y and z directions, respectively. 

Another nxamerical stability criterion for the free-surface 
model is the Courant-S^iedricks -Lewy (CFL) condition, Roache (1972)^ 
and Richtn^yer. and Morton (1976), which is based upon external 
gravity waves (or so-called surface gravity waves) and is expressed 
as follows : 

At < — or At < ^ (3.41) 

Vp v/P 
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whichever is smaller, is defined as the phase vel- 
ocity or the celerity of these external gravity waves , 

3,7,6 Nximerical Modelirtia: Approximations 

3 , 7 , 6 , 1 Adi^atic Condtion 5.n FarT.Field' Version 


Following Roache C1972) , the lateral temperature boundary 
conditions is specified as: 

T = T ,1 on 0 = constant planes (3.42) 

w w +L 


This condition is simply set after the interior point cal- 
culation for Vi is performed , Although 


11 = .iT,13_n-0yis 

3a 3a^H3a H3a-' 


(3.43) 


and , 

^=3JD + 3_T.l^-ayi>. (3.44) 

3y 3g 3a '' H 36 H 3 3 

9 T 

for a shallow body of water like the Biscayne Bay, is 

quite small in comparison with the horizontal temperature var- 
iations, and, therefore, is neglected by using on 

constant planes. 


3 , 7 , 6 . 2 Velocity Gradient at Inlet, Far-Field Version 

The velocity gradient (or as the case might be) 

has been approximated as follows (for Biscayne Bay) 

16 

V (I,J,K) = 2 V (I,J^,K)/10 (3.45) 

inlet interior 
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where corresponds to the value of the J index at the interior 
points adjacent to the inlet. This approximation has been used, 
since only one value of tIq (t) is known at the inlet for the 
Biscayne Bay, However, this would not be essential if rioCt) was 

known at all points along the inlet’ 

3,7,7 Higher Order Terms 

Higher order terms resulting from the transformation of 
the horizontal diffusion terms Ci.e, second order derivatives 
in a and B ), from the Cartesian coordinate system (x,y,z) 
to the vertically stretched coordinate system Ca,3,c;) have 
been neglected. This has been done since the magnitude of the 
vertical diffusion terms are several orders of magnitude 
larger than the horizontal diffusion terms, Sengupta and Lick 
(1974) , Appendix A of this volume presents the details of this 
transformation , 

3 , 8 Sample Results 

In this section sample results using the free-surface 
model for near-!-field and far-^field applications will be pre- 
sented. 

3.8.1 Far-Field’ Biscayne Bay (closed basin with ocean efflux) 

Ignoring the Cutler Ridge site thermal discharge, the 
free-surface far-field model was applied to Biscayne Bay to 
investigate the general circulation, natural temperature dis- 
tribution and the surface height behavior. Lee and Sengupta (1977) . 

Preliminary cases and governing physical factors were 
first studied. Table 3-1 gives the various cases run, Model 
execution including all the physical factors, wind, current, 
tide, bottom topography for the April 15, 1975 data base was 
the final case run for calibration and verification. 
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Fig^ 3,5 shows the model predicted surface velocity dis- 
tribution at 0900 EST on April 15, 1975, This case is for 
incoming tide with the tidal current velocity specified 
sinusoidally with respect to time at the oceanr«bay interface, 

The -effect of the wind is clearly indicated to be in the dir- 
ection of the wind in the northern closed part of the bay, 

Fig, 3,6 shows the model predicted surface velocity distribution 
at 1300 EST, This case is for outgoing tide and the south 
wind effect is clearly exhibited, In Fig, 3,7 the velocity 
distribution at a depth of 1 meter is shown where it can be 
seen that the tide effect dominates the effect of the wind, 

Fig, 3,8 shows the contours of constant surface height 
at llOOEST as predicted by the model , and Fig , 3,9 shows the 
contours of constant surface height at 1400EST, Now Fig, 3.10 
and Fig, 3,11 show the corresponding surface velocity distribu- 
tion at 1100 EST and 1400 EST, respectively. It can be seen 
that comparison of Fig, 3,8 with Fig 3.10. and comparison 
of Fig, 3,9 with Fig, 3,11 indicate the relationship between 
the lines of constant surface height and the velocity field, 

. s 

Fig, 3,12 shows two synoptic model isotherm plots ya 
IR-data for April 15, 1975, The agreement is good, Fig, 3,13 
shows the surface height versus time at two tide gaging 
stations, observed vs calculated. The agreement is relatively 
close, Fig, 3,14 shows surface height versus x-direction 
Cl-direction in grid system) for varying time during the tidal 
cycle along the transect J=7 , It can be seen that the surface 
moves according to the stage of the tidal cycle, as 
would be expected, 

3,8,2 Combined Near«?Field and Far<^Field by Using Hor- 
izontal Stretching for Hutchinson Island (open basin) 

Fig, 3,15 shows the general location of the Hutchinson 
Island Power Plant site. The condense cooling water for the 
power plant is provided by the intake and discharge pipes which 
circulate the water from the Atlantic Ocean with canals to the 
plant, Fig, 3.16 shows the details of the submerged discharge 
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pipe geometry. 

Fig. 3. 17 shows the horizontal grid system without hori- 
zontal stretching (I“20, J«20) . Fig. 3. 18 shows the horizontal 
grid system with SINK stretching. Fig. 3. 19 shows a distorted 
vertical section with a- stretching. 

The model was run using the June 2, 1976 data base. On 
June 2, 1976 the thermal discharge had an exit velocity of 
280 cm/ sec (9,1 ft/sec) at each end of the Y-type discharging 
pipe. The discharge temperature was about 35°C. 

Fig. 3. 20 shows the surface height variation along the 
1*8 transect (one grid point before the discharge points from 
south to north) . The storface heights are negative around the 
area of discharge owing to viscous entrainment. Fig. 3. 21 
shows the surface velocities with the conditions for June 2, 1976. 
The imposed northerly current prevails away from the discharge. 
Fig. 3.22 shows the horizontal velocity on the plane of the 
discharge. The velocity field is dominated by the discharge 
conditions because the inertia of the jet is the important 
driving mechanism near the discharge point. The velocities 
decreased away from the discharge owing to entrainment as is 
expected. Fig. 3. 23 shows velocity distribution in the J*10 
transect which is perpendicular to the shoreline. A small 
vortex can be seen near the discharge region owing to entrain- 
ment just west of the discharge. 

Fig. 3. 24 shows the comparison of model results and IR 
data for June 2, 1976. Relatively good agreement is observed. 

The model was verified for May 17, 1977. The free surface 
model was run for one hoxir with May 17, 1977 data base as an 
input. Fig. 3. 25 shows the surface isotherms comparison of model 
results and IR data (1113-1118 EDT) . The isothems of 25.4°C 
and 25.9°C from model results cover a larger area than IR data. 
Generally, the results predicted by the model are in agreement 
with the IR-Data. 

The results presented are taken from Tsai (1977) . This 
model is extremely sensitive to parameters of the problem, 
time step and boundary condition. Further verification is 
ongoing. Where high resolution at discharge point is not 
necessary the horizontally stretched free-surface model is 
recommended. 
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TABLE 3-1 LIST OF COMPUTER CASES RUN FOR 
FREE -SURFACE FAR FIELD - BISCAYNE BAY APPLICATION 
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**Only results for this computer case are being presented at this time 























TABLE 3-2 LIST OF CASES OF THE FREE SURFACE MODEL APPLIED ' 
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Fig. 3-5 Surface Velocity Distribution for Biscayne Bay 

at 9:00 AM, .April 15, 1975 (by Free Surface Model) 
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Surface Velocity Distribution for Biscayne Bay at 
1:00 p.m., April 15, 1975 (by Free Surface Model) 
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Fig. 3-10 Surface Velocity Distribution for Biscayie Bay 

at 13:00 am, April 15, 1975 (by Free Surf ice Model) 
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Fig. 3-11 Surface Velocity Distribution for 

, Biscayne Bay at 2:00 p.m,, April 15/ 1975 
(by Free Surface Model) 
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Fig. 3-13 Surface Height vs. Time, at 

Miami Beach and Coconut Grove, 
Observed versus Calculated 
(Free-Surface Model) 
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il Grid Point System Without Stretching 
Surface Near Field Model Applied to 
n Island Site 
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Fig, 3-18 Hori, 2 ontal Grid Point System witl* SINH Stretching 
For Free Surface Near Field Model Applied to 
Hutchinson Island Site 
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Fig. 3-20 


Surface wave height along 1=8 from south to 
north at Hutchinson Island Site (Free Surface 
Model) 
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Pig.. 3-21 Surface velocity distribution with current, 
wind and bottom topography at Hutchinson 
Island Site (Free Surface Model) 
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Fig, -3-22 Horizontal Velocity Distribution near the 
ocean bottom around the discharge pipe at 
Hutchinson Island Site (Free Surface Model) 


146 


THERMAL POLLUTION LAB 
UNIVERSITY OF MIAMI 


9 Z 

O O 09 

0 ) o> q 

SI 

fi a 

q 8 8 q 
N q *H 


0^0 

u 

00 • rH lO a 

evi ^ w 

a > 

1 1 

1 1 

>k 

«p 

>% 

•r4 

A 

q 

Oi 

O 

a 

H 

u 

P 

bO 

> 

& 

q 

p 

bo 

H 

p 


cil 

a a 


q 0 

q 

u ^ 

m q 



3 0 


on 


i 


f 


A 


4 

1 

k 

f 


A 


f 


A 



CO 

q 

O -rl 
0) -P 
C 

•H \n 

^ Ol 

o • 

4J in 


u -P 
o OJ 




I 

b 


g 


M 

q 

01 

01 

(d 

X 

q 


q 


4J Id 

o o 
q n 
q 

H 
•P (d 
Id o 


g 4J • 

o 

•H q H 
4J > q 
0 ^ 
rQ rC 0 

•H 4J g 
P *H 

q 

w o 

•H q Id 

-P m 
•H p 
>iw :3 
•p w 

•H ^ 

o c q 

o q q 

H rH p 

q W fl4 

> 


on 

CV| 

I 

on 


01 

•H 

U4 



THERMAL POLLUTION LAB 
UNIVERSITY OF MIAMI 


/I 

/ 

/ 

> 

y 

j 

/ 

j 

i 

i 

/ 

y 

s 

y 

y 

y 

y 

y 

y 

y 

y 

y 

y 

y 

y 

y 

y 

y 

y 

y 

y 


Pig. 



Discharge Temp : 
Air Temp : 
Ocean Temp : 
Current : 
Wind : 

Bottom Topography: 


35°C 

29°C^ 

25.5°G- 
25 cm/sec N. 
4.47 m/sec 
(10 mph) S.E. 
Varied 



Result 


0 250 

U » J 

Distance (meter) 


3-24 Comparison of model results and 
afternoon IR data at Hutchinson 
Island Site for June 2, 1976 
(Free Surface Model) 
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Fig. 3-25 Comparison of model results and IR 

data (1113 - 1118 EDT) at Hutchinson 
Island Site on May 17, 1977. (Free 
Surface Model) 



IV, EVALUATION OF MODELS 


Verification of existing models has been in general 
quite unsystematic, In order to provide a basis for user 
confidence it is essential that accuracy and versatility 
of models be established by repeated verification at diverse 
sites . While the verification procedure for the present 
efforts is by no means completely satisfactory, significant 
improvements have been made in the thoroughness of verifi- 
cation, A significant factor contributing to this improve- 
ment hao been the integration of remote sensing and insitu 
data acquisition program with the model development effort. 
Accurate specification of initial conditions has been greatly 
enhanced by synoptic data bases for IR measurements . This 
have been a serious drawback of other efforts in thermal 
pollution model development to date. 

Comparison of performance of different models is diffi- 
cult to make. There are two reasons for this difficulty. 

1. Systematic synoptic data bases which can be used 
as standards do not exist, 

2 , Performance of a given model is dependent on the 
validity of assumptions and approximations for a given site. 
Thus performance is often site specific, 

An attempt at comparative evaluation was made by Dunn 
et al (1975) . Their extensive review has been documented 
in a two volume report. Some models were compared using 
standard data bases for che Point Beach Plant. However, 
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confusion regarding initial and boundary conditions were 
still present. Since the models developed by the University 
of Miami team were calibrated and verified at locations where 
no other model has been used, detailed comparative evaluation 
is not possible. However, in order to present the performance 
of the models developed, in perspective, results presented for 
various models by Dunn et al (1975) are discussed briefly. 

Table IV- 1 shows a summary prepared by selecting the 
most commonly used models. These encompass phenomenological, 
integral and nxomerical models. Details of these models with 
critical comments are presented by Dunn et al (1975). Fig. 

4-1 to 4-3 shows comparisons of results from models by Waldrop 
and Farmer with data at John Sevier Plant and Point Beach 
Plant. Surface isotherm predictions at the first site are 
relatively good. Errotsof approximately 2°C are present in 
vertical temperature distributions . Comparison of centerline 
trajectory is poor. Centerline temperatures show large errors 
after 200 meters. Area under isotherm, predictions fall ap- 
proximately 7 time below measured values . No velocity com- 
parisons were made, 

The Stolzenbach-Harleman integral model was tested for 
results at different tide stages. Fig. 4-4 to 4-6 shows sur- 
face isotherm predictions. The model consistently under pre- 
dicts isotherm areas . The comparisons being poorest at low 
tide indicating that bottom topography effects are not ade- 
quately modelled. 

Prichard's phenomenological model was compared to the same 
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data base as the Stolzenbach-Harleman model. Figs 4-7 to 
4-9 shows comparisons at high tide, mid-tide and low tide. 

The isotherm areas are under-predicted, with error being 
maximum at low tide. The results are qualitatively better 
than the Stolzenbach-Harleman model. However, bottom topo- 
graphy effects make the model non-usable in practical 
applications . 

Comparisons of results using Prych's phenomenological 
models are shown in Figs 4-10 and 4-11 for applications at 
Point Beach Power Plant and Waukegan Power Plant. Relatively 
good agreement is observed for the Point Beach case except 
very close to the discharge point. The comparison for the 
Waukegan Power Plant is significantly worse. 

One of the models that can be used relatively easily 
is the one presented by Shirazi and Davis (1974) . They 
present nomograms and sample problems in a two volume work- 
book. Figs 4-13 and 4-14 show comparisons of predicted 
values and mean data from a number of sources the agreement 
is good. However, this model is quite unsuitable for basins 
where the infinite depth assumption is not valid. 

The ntomerical model of Till (1974) has been compared 
to field data obtained at Phillip Sporn Power Plant. Figs 
4-15 to 4-17 show isotherms in vertical sections . Near the 
discharge an error of about 2°C is observed. Comparisons 
become better with distance from discharge. 

The Paul and Lick (1974) model is very similar to the 
rigid-lid model developed by the NASA-KSC, University of 
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Miami effort. The comparison with Point Beach field data 
are shown in Fig. 4-18 and 4-19, Underprediction of areas 
under given isotherms is observed. This could be owing to 
errors in choice of diffusion coefficients. Difference in 
plume shape may have been caused by inadequate information 
regarding ambient currents. 

The comparisons for the present models have been pre- 
sented before. The following summary statements can be 
made, 

a) . The comparisons of rigid-lid near field model for 
Cutler Ridge plxme is in good agreement with IR data as 
shown in Fig 2-11, The centerline temperatures are espe- 
cially well predicted as shown in Fig 2-12. 

b) . The predictions of rigid-lid model for Lake Belews , 
mixing pond is in agreement with IR data to within 0,2°C 

as shown in Fig 2-30, The main lake predictions have shown 
errors of upto 3°G at narrow cross sections owing to lack 
of spatial resolution as well as uncertainty in data regar- 
ding the thermachine location, 

c) , Comparisons of free surface model results with 
field data at Hutchinson Island show good agreement both 
for plume shape and temperature as shown in Fig. 3-24. 

d) . The far-field rigid lid model applied to Biscayne 
Bay shows surface isotherm predictions to be within 1°C 

of corrected IR data, as shown in Fig 2-22, 

e) . The far-field free-surface model predictions for 
Biscayne Bay show agreement to within 1°C. 
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It is imperative at this stage to note that Little or 
no velocity verifications exist for any of the models devel- 
oped. This is true for all the models presented by Dunn 
et al (1975). In the present study some limited float mea- 
surements of surface velocities in Lake Belews were obtained. 
Qualitative agreement with model results were observed as 
reported by Mathavan (1977) , For complete verification of 
models, velocity verification is essential. However, until 
field measvirement equipment that can accurately measure 
velocities in the range of 0 - 10 cm/ sec is developed, such 
verifications cannot be made. 
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Fig . 4-3 Comparison of Waldrop-Farmer Model Predictions to Two Sets of Temperature Data 
Taken a: Point Beach (Unit 1) on September 1. 1971. 
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FIELD DATA 


Fig r 4-15 Isotherms from Field Data and Calculated Values for Philip Sporn 
Power Plant at Position C, Approximately 1000 ft Downstream 
from Discharge Point. (Adapted from Till, *) , 


CALCULATED VALUES 


Fig. 4- 1:6 , Isotherms* from Field Data and Calculated Values frr Philip Sporn 
Power Plant at Fositum D, Approximately 3000 ft Downstream 
from Discharge Point. (Adapted from Till. ) 
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Fig . 4- 17 • otherms from Field Data and Calculated Values for Fliilip Sporn 

-owti Plant at Position E, Approximately 4400 ft Downstream 
from Discharge Point. (Adapted from Till.^) 
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t 

Fig . 4rl8 Comparison of Surface Isotherms from Measurements and Predictions 
of PauI^Iick Model: Point Beach Power Plant, May 18, 1972, 



Fig . 4-19 Comparison of Surface Isotherms from Measurements and Predictions 
of Paul-*Uck Model: Point Beach Power Plant, May 20, 1971. 
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V, CONCLUDING COMMENTS 


From the experience of existing modelling efforts it 
can be concluded that numerical models are the only viable 
tools that can successfully incorporate the physical details 
of discharge, receiving basin and atmospheric conditions. 
While phenoinanblDgical and integral models are easier to use 
they contain generic deficiencies that cannot be solved by 
continued calibration and refinement. Therefore, the search 
for models which can be used with confidence for regulatory 
purposes should be directed to numerical models only. 

The model package developed by the University of Miami- 
NASA-KSC efforts show the promise of providing a reasonably 
general model package. Further calibration and verification 
of these models should inspire user confidence. 

Some comments regarding research efforts that are needed 
for the development of generally accepted models is appro-' 
priate to conclude this volume, 

1, Formulae for eddy transport coefficients need to 
be developed and verified, 

2, Surface heat exhcange coefficients and radiative 
transport of heat into aquatic domains need to be better 
understood, 

3, Reliable anemometers sensitive at low velocity 
ranges need to be developed. 

4, The surface "skin” temperature profile must be 
understood in terms of meteorological and surface turbulence 
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conditions, before extensive used of IR data can be made 
without ground truth measurements, 

5. Most numerical models are relatively expensive to 
use owing to computer time costs. However, the cost of 
computer time is minor compared to overall cost of environ- 
mental impact statements. The situation can further improve 
through development of more sophisticated numerical methods 
as faster computers, 
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APPENDIX A 

DERIVATION OP STRETCHED EQUATIONS 
[I] Vertical Stretching Equation 

To incorporate with both free surface and realistic 
variable bottom topography in the three-dimensional basin 
model is quite difficult in computer progrzunming . A mathe- 
matical transformation of vertical coordinate is needed to 
convert the depth of basin to a constant non-dimensional 
depth one. The new coordinate system is transformed from 
X, y, z, to a, 0, a, where the transformation relationships 
are: 

^ _ Z (g,8,z,t) ^ z + n (q#B#t) t , 

H H (a,0) + n (<x,6,t) 

ax ay vA z) 

If - If - H - H- “ 

H - 1 <A-4) 

= i, (A—K\ 

az . az H * ' 

a g _ ^ 3n _g an 

ax “ H 3a H 3o ’ 

3 0 _ i 3 n _ 2. iA—i\ 

ay ~ H as H 30 ' 

By using the above relationships the first derivatives 

can be written as 

+ i£_l£ 1£ 

3x ~ 3a ax 30 ax 3cr 3x 3ct 3a 3x“3a *■ 3a 

rian. lisill (a a\ 

^ H 3a “ H 3a ^ " 3o H 3o 3a ” H 3a 3a ^ 
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(A- 9) 


lU.m + ilalU £M1£ 

3y 36 H 30 aa " H 36 3o 


iZ « 1 

3z H 3o 


(A-10) 


Where "F” is. the appropriate dependent variable/ the 
second derivatives can be written as: 


3^F , 3^F _2^^aa.2aa. 3^F ,30v 2 32p 

3x2 “ ff 3o da ax 3X 3aaa 'ax^ a^ 

i 1£ 

H 3ff 3a2 H 3 ct 3a 2 * 


(A- 11) 


3^F ^ a2p 2 ^ 3H 3a . 2 30 3^p /i2.x ^ 

3y^ WS" ” H 3< 36 3y 3y 363cr ^3y' 3a2 

4. i 1£ iio. £ IF 

H 30 302 " H 3<T 302 

a2p 1 a2p 

3i? * if2 W? 


Define : 


U 


dx _ do 

at " at 


V 


dv de 

at ' at 




da 

■ at 


(A- 12) 
(A- 13) 

(A- 14) 
(A- 15) 
(A- 16) 


The vertical velocity w is related to the fi by the 


expression 

dz Ty it 1 \ dr) . dh 
w - - Hn + (cf - 1) dt at 


(A- 17) 


Continuity Equation : 

The continuity equation for three dimensional incom- 
pressible fluid flow in Cartesian coordinate may be written 


as : 




H 

« ay 


+ H = 0 (A-18) 

d Z 


178 




' w 


Is 


Now 


„ au „ „ 3u . In 5u ^ 3H au 
“ax*®^+aa37-®mi7- 


„ av _ „ av , an av ^ an av 
“ay “ a$ ■’■ ae ao “ ® a? To’ 


(A-19) 
(A- 20) 


n 3'' 
“ ■?* 


aw ^ 3w _ a /dZ. _ a / d(Z - n ) 1 
“ -sr ■ ~ 'Tnr) ■ ' dt ^ 


“ 7J “ 3o 3o '3E 


a .dZ _ 3 /d(oH) . dnv 

3o '3E " 3t^ 3a ' dt 3t' 

a /i, do , „ dH _ dn > 

^ at + ® 3 E at^ 

“ a 7 <dt^ at ° 37 <at^ ^ <at> 


Hia + /3H + u^ + v^+niiii + o-^(^ 

“ 3o ^ '3t ^ “ 3a ^ ^ 3B ^ ” ^o’ ^ ® 3o 'at 


+ U 


3H 

3a 


X „ „ 3 H. 3 /In X « in X V in. + n 


„ an . 3H . „ 3H . „ ^ ^ au an . „ 3v ^ 

“ 3o at 3o ^ ae ® ao ao ao ae 


au an av an 

” 3o ao 3a 30* 


(A-21) 

Substituting (A-19 ) , (A-20 ) , and (A-21) into (A-18) we get 
the continuity equation in a 0 o coordinate is : 

0 (A-22) 


iS X 3 (Hu) . 3(HV) X H ^ 
at 3a 30 ao 


Momentum Equation and Energy Equation 

The inertia terms in the momentum or energy equations 
in the xyz Cartesian coordinate can be written as ’ 


H H - I - If ^ H- 


(A-23) 


where P is the appropriate dependent variable u, v or T for 
the momentum and energy equation respectively. By using the 
relationships in (A-8) , (A-9) , (A-10) , (A-17) and added F 
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1 


times the continuity equation [P (|^ + ^ + H |^) ] , 


3a 


3$ 


3o' 


then the equation (A>23) becomes: 


ilHH ^ 3JHFU) ^ i l mrjO H liPQl (A.24) 

Ob ooi 33 vO 

The pressure terms can be written as: 

F 9x F 3a ^ 3o 3o " pH To 3o 


1 3P 
^ 3a 


1 3P 


3H 


Jii QF A OK- / 0 n \ 

" ^ 3a 3a “ 


_ 1 3P 1 3(p q Z) ,_3H Bha 
p ^ ^ ^ 

^ 3a ” H 3o ' 3o " 3a' 

-iff- s(» If - If) 


1 3P _ 1 3P 3H 3n» 

pT 3y ? 30 " 30 “ 30^ 


(A-26) 


By using the relationships of second derivatives, the 
diffusion terms can be written as: 

3^F ^ 2 3o aja ^ 2 30 3^F 2 a2p i 

3x2 3a2 " H 3X 3a 3a 3x 3a3o '3x^ 3o2 H 


il. 3^n _ 0 3F 3^H 
30 3x2 “ H 30 3o2 

1 32p ^ ip. 1 3H ^ 

“ H 3^ 3o 3a' “ H 3« 3a 


a 21 3^H 

H 30 


2 ^ 3H ^ . 2 30. 3^F . ,30» ^ 3^F . 1 3P 
H 3x 3o 30 3x 3a30 'Sx' H 30 3^ 


[ 37 (“Ifn + High order terms, 


(A-27) 


using the same procedure as above: 

0 - ff is? H’) ’'iSh 


(A-28) 


and the vertical diffusion term is: 
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Naglectlng the high order terms in Equations (A-27) and 
(A~28) , then 

1^ " ff <A-30) 

||t ■ ff tii |f)l <A-31) 

By using the relationships in Equations (A-24) , (A-25 ) , 
(A-26 ) , (A-29) , (A- 30) 2 md (A‘*31) / then the momentum and 
energy equations become: 
u~momentum: 

3 (Hu) . 3(HUU) ^ 3(HUV) . 3(uft) 

^ 3a * H “ ”37~" 

- B t- i (||) + 9<o I? - I5-) + IV] + K„ (H |a>] 


Us as 


v^momentum: 


3a 


3a' 


'H ^3a 


3a' 


>1 + F 'hI? '■> Bv l?)l 


3HV . 3(Huv) , 3(Hw) , „ 3(vfl) 

"Tt — 33 — “ “T7“ 


B t- F <H' 


B*® 86 “ 86*' Uo *B a„)) + Bgtjj (H ) 

tsf? <» Bvl?)] <A-33) 

Energy Equation : 

3(HT) . 3(HuT) . 3(HVT) . „ 3(£2T) _ „ r 3 /« 8T» i 

Tt'- + 53~ * 86 * B -JS “h 'B 3„) J 

‘“H'l 'gl? ‘O BvH>l 

[II] Horizontal Stretching Equation 

Before writing the equations in finite difference form, 


the horizontal stretching transformation which provided a 


I ' 


r' 

\ 


niortt •fficlant u«« of grid points within a large domain 
was applied In a and 0 direction. The hyperbolic sine 
stretching equation was used by letting: 


a - a + Slnh (Cj (X - d) ] (A-35) 

0 ■ b + Cj Slnh tC4 (Y - e)] 1 (A-36) 


where a,0 Is a real coordinate, X and Y are stretched 
coordinate, a and b are the distance at which the minimum 
step size Is desired. ^2 ^3' ^4' constants 

to be detexnoined by the Imposed conditions. 

The differential derivatives transformation relation- 
ships for the new X Y a coordinate from a 0 0 coordinate are 
written as: 


3P . 3X ap _ 3P 
^ ^ 33f “ * 33? 


(A- 37) 


3P _ 3Y 3P _ 3P 
30 30 3Y ^ 3Y 


3^P ^ _3 
3a2 9a '3o' 


(X*)^ 



X" 


1£ 

3X 


3^P _3 . 3Pv 

302 90 '90^ 


(Y')2 


9^P 


9P 

•5Y 


(A-38) 


(A- 39) 


(A- 40) 


where 


X’ 


12 . Y* » 

3o ' ^ 30 


X” 


3^X 

3o7 


Y« 


a^Y 

307 


By using the above relationships, then tne set of 
equations in transformed X Y a coordinate are: 
Continuity Equation ! 


3H 3 (Hu) 3(Hv) , „ 30 _ „ 


(A-41) 


-1 
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Mom»ntum Equaticais : 


or 


• Page is 

QUAuir 


u-moncntumi 

9 ^ + X. »^Baai. + Y- 91 |§ 2 L + H . H (-S’ (||) 

* » <» II - ■+ *''J + ■'b '«■>* I H 

+ H (x')»|^ + HX*|5] +Kg ((y)9|?|§ 

+ H(V)9 + H V |§J + 1 [g|^ (p K, |a)].. (A-42) 

v-nomentum: 

3 (Hv) . 3 (Huv) y, 3 (Hw) .. 3(vfl) 

3t * 33C * “T5f ”^“30 

■ “I- V ‘ 1 ^* « ’'■ <“ H - tJ’ ■ 

* Kg((X*)9 HIJ+ B(X')9 |^ + HX* IJl 

* #W + 

*? IhI^ <'■ ^I 5 >i <^-«> 

Energy Equation ; 

3(HT) . 3 (HuT) . 3 (HvT) . „ 3(nT) 

"It * 3X ^ 3Y “ 3a 

■ HH 1 ^* “ =‘" H> i 

*»H If?* “’'"Wl 

*1 IhI? '» «vH'l 
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APPENDIX B 

Hyperbolic Sine Horizontal Stretching System 

It Is desirable to obtain a more detailed description 
of the flow near the discharge point while large grid size 
may be used In the further points to prevent the unnecess- 
ary computation time. A horizontal stretching has been In- 
vestigated and Incorporated to allow for high resolution 
near the discharge with gradually Increasing grid sizes aw- 
ay from the discharge. The tangent equation which was used 
by Waldrop and Farmer (1974) appeared to have the most des- 
irable characteristics. However, it has been noted that 
when the number of grid points is small and the domain is 
large, the tangent equation produces relatively little str- 
etching until close to a boundary, then jumps to the bound- 
ary in a relatively large steps. This behavior appears un- 
desiarble and further effects to find a suitable alternative 
equation. The hyperbolic sine has been investigated by C. V. 
Carter (1976) . It appears to have somewhat better charact- 
eristics than the tangent equation. 

(I) The Hyperbolic Sine Stretching Equation : 

The hyperbolic sine stretching equations used for both 
transverse and lateral directions are 


X - a + Cj^ Sinh{C 2 (X-d)} (B-1) 

y - b + C 2 Sinh{C^(Y-e)} (B-2) 


where x and y are the real horizontal coordinate; a nad b 
are the distance at which the minimum step size is dejired; 


^1' ^2* ^3* ^4' con*tant8 to b« daterminad 

by tha imposad conditions; X and Y ara tha atratchad horiz- 
coordinata. 

Tha variables X and Y ara computed as 


X - (I-l) AXq (B-3) 

Y ■ (J-1) Ay^ (B'-4) 


where I is tha grid point number on x-axis; J is the grid 
point number on y-axis; Ax^ is the minimum desired step size 
in x-direction; Ay^ is tha minimum desired step size in 
y-direction. 

In order to determine the constants, we impose the 
following conditions : 

(a) When x « 0 then X ■ 0, y ■ 0 then Y ■ 0. Form equati- 
ons (B-1) and (B-2) , this condition will be satisfied if 


1 1 

d * Sinh”^ (-^) (B-5) 

^2 ^ 

e “ — — Sinh ^ ( ^ ) (B-6) 

C 4 ^3 

(b) When it reaches to the boundary, the equations (B-1) 
and (B-2) have the form as 

X|^ * a + Cj^Sinh{C 2 CNjj,-1)AXq - d } (B-") 

yb » b + C 3 Sinh{C 4 (Ny-l)Ay^ - e } (B-8) 


where x- and yb are the x and y reach to the boundary resp- 
ectively; N„ is the total number of grid points on the x-axis ; 
Ny is the total number of grid points on the y-axis ; then the 

c»-3 
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valu«8 of X and Y at the boundary are (Njj-I)AXq and (Ny-l)AyQ 
respectively. 

(c) Both (a) and (b) above are the required boundary con- 
ditions. The imposed conditions at x « a and y "■ b which can 
be anywhere in the domain. When x ** a, it is required that 
the step size be minimm in x-axis. When y ■ b, it is req- 
uired that the step size be minimum in y-axis. That is 


Ax “ AXq when x»a (B-9) 

Ay - Ay^ when y ■ b (B-10) 

one can write 

Ax - AX (B-11) 

Ay - AY (B-12) 

but AX « Ax„ ; AY - Ay^ , therefore 

0 •' o 

** ""3X" ’ * (B-13) 

Ay - Ay^ (B-14) 

From (B-9), (B-10), (B-13), (B-14), we can find that 

» 1 when x«a (B-15) 

*» 1 when y ** b (B-16) 


Differentiating equation (B-1) and (B-2) and setting the 
result equal to 1, we can find that 
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^2 ** (B-17) 

Ca ■ — ^ (B-18) 

3 

Substituting equation (B-17) into (B-1) and (B-18) into 
(B-2) and the result are 


X - a + Sinh( ) (B-19) 

y - b + C3 Sinh( ) (B-20) 

3 

By substituing equation (B-17) into (B-5) and (B-18) into 
(B-6) , we can get 


d 

e 


C3_Sinh"^(-gj) 

C3Sinh"^(-^) 

3 


(B-21) 

(B-22) 


At the boundary, the equations (B-19) and (B-20) becomes 

(N -l)Ax - d 

Xb - a + Cl Sinh {— ^ } (B-23) 

(N -l)Ay - e 

y^ *• b + C3 Sinh{ — ^ } (B-24) 

3 

There is only an unknown in the equations (B-21) 
and (B-23) , an unknown C3 in the equations (B-22) and (B-24) . 
The iteration method was used for solving the equations to 
obtain and C3. The d and e are obtained by substituting 
and C3 into equations (B-21) and (B-22) respectively. 

The derivatives -gl- j 5 -3^ ; are needed 

for the governing equations when horizontal stretching equ- 
ation are used. The derivetives are obtained from differen- 
equation (B-19) and (B-20) and that are: 
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dX _ 

usr " 


d^X 

Tx2" 


dY 

w 


d^Y 
dy ^ 


X' - 1 (B-25) 

Coah 

Sinh 

. X” - ~ 4 ^ (B-26) 

^ Co«h»(^) 

Y' « 1 (B-27) 

Cosh 

*^3 

Sinh 

, Y" - - 4 — (B-28) 

"^3 Cosh 3 (L^) 

^3 


(II) Numerical Results : 

The following characteristics were chosen as an init- 
ial domain for the Hutchinson Island Site: 

Y 


A 

Y - Boundary 

^ 

1 

1 

1 

1 

1*^ 
' 1 

Discharge ja 

Points 2 

! ^ 
1 y 

1 ^ 

t 

1 

1 


i<-.a 
' 1 
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TABLB 3-1. Samp] 

Le Of Numerical Results 


X-Axis 

Y-Axis 

Total number of points 

20 

20 

Boundary distance 

Xjj m 238000cm 

y^ * 200000cni 

Discharge pipe outlet 
location 

a ■ 38000 cm 

b ■ 100000 cm 

Minimum desired step 
size 1 

X — 1500 cm 

Y » 1500 cm 

Constant 

m 3696.57 

Cj * 3530.03 

Constant 

d « 11184.66 

e * 14251.88 


Pig,B-l and Pig.B-2 are the comparison of hyperbolic 
sine stretching equation and tangent stretching equation. 
These graphs show that for the specified conditions, the 
tangent equation produces relatively little stretching until 
quite close to the boundaries. The sinh equation is somewhat 
better in this respect and stretching occurs more gradually 
as X and y vary between their limits . 
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